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Abstract 


The contents of this report covers: (i) development of 
optimal geometry for crowned helical gears; (ii) method for their 
generation; (iii) tooth contact analysis (TCA) computer programs 
for the analysis of meshing and bearing contact of the crowned 
helical gears and (iv) modelling and simulation of gear shaft 
deflection . 

The developed method for synthesis is used for determination 
of optimal geometry for crowned helical pinion surface and is 
directed to localize the bearing contact and guarantee the 
favorable shape and low level of the transmission errors. 

Two new methods for generation of the crowned helical pinion 
surface have been proposed. One is based on application of the 
tool with a surface of revolution that slightly deviates from a 
regular cone surface. The tool can be used as a grinding wheel 
or as a shaver. Other is based on crowning pinion tooth surface 
with predesigned transmission errors. The pinion tooth surface 
can be generated by a computer controlled automatic grinding 
machine . 

The TCA program simulates the meshing and bearing contact of 
the misaligned gears. The transmission errors are also 
determined . 

The gear shaft deformation has been modelled and 
investigated. It has been found the deflection of gear shafts 
has the same effects as those of gear misalignment. 
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SUMMARY 


The topology of several types of modified surfaces for 
helical gers are proposed. The modified surfaces allow 
absorption of linear or almost linear function of transmission 
errors caused by gear misalignment and deflection of shaft. 
These surfaces result in the improved contact of gear tooth 
surfaces. The principles and corresponding programs for computer 
aided simulation of meshing and contact of gears have been 
developed. The results of this investigation are illustrated 
with numerical examples. 

1. INTRODUCTION 

Traditional methods for generation of involute helical gears 
with parallel axes provide developed ruled tooth surfaces for the 
gear teeth (Fig. 1.1). The tooth surfaces contact each other at 
every instant along a line, L, that is the tangent to the helix 
on the base cylinder. The surface normals along L do not change 
their orientation. The disadvantage of regular helical gears is 
that they are very sensitive to misalignments such as the 
crossing or intersection of gear axes. The misaligned gears 
transform rotation with a linear function of transmission errors 
(a main source of noise) and the bearing contact is shifted to 
the edge of the teeth. The frequency of transmission errors 
coincides with the frequency of tooth meshing. The actual 
contact ratio (the average number of teeth being in mesh at every 
instant) is close to one and is far from the expected value. 
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These are the reasons why we have to reconsider the 
canonical ideas on involute helical gears and modify their tooth 
surfaces. Crowning of the gear surfaces is needed to negate the 
effects of transmission errors and the shift of contact between 
the gear tooth surfaces. Deviations of screw involute gear tooth 
surfaces to provide a new topology that can reduce the gear 
sensitivity to misalignment will be developed. Theoretically, 
the modified tooth surfaces will be in contact at every instant 
at a point instead of a line. Actually, due to the load applied 
between meshing teeth, the contact will be spread over an 
elliptical area whose dimensions may be controlled. Methods for 
gear tooth surface generation that provide the desirable surface 
deviation are proposed. For economical reasons only the pinion 
tooth surface is modified while the gear surface is kept as a 
regular screw involute surface. 

2. BASIC CONCEPTS AND CONSIDERATIONS 
2.1 Simulation of Meshing 

The investigation of influence of gear misalignment requires 
a numerical solution for the simulation of meshing and contact of 
gear tooth surfaces. The basic ideas of this method (Litvin, 
1968) are as follows: 

(1) The meshing of gear tooth surfaces is considered in a 
fixed coordinate system, S f . Usually, the generated gear tooth 
surfaces may be represented in a three parametric form with an 
additional relation between these parameters - Gaussian 
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coordinates. Such a parametric form is the result of 
representation of a gear tooth surface as an envelope of the 
family of the tool surface (the generating surface) and two from 
the three Gaussian coordinates are inherited from the tool 
surface . 

The continuous tangency of gear tooth surfaces is 
represented by the following equations 


r ( u l f0 i' 1 J'i , 4 , i) £ ^ u 2'®2 , ^2 ,< ^2^ 


( 2 . 1 . 1 ) 


( 1 ) 


( 2 ) 

— ' • 


n (u 1# e 1 ,* 1 ,* 1 ) = n v (u 2 ,e 2 ,^ 2 ,4> 2 ) , 


(1) 


(2) 

n 

— 

n 


( 2 . 1 . 2 ) 


f 6 (u 1 ,0 1 ,^ 1 ) = 0 (2.1.3) 

f 7 (u 2 ,0 2 ,^ 2 ) = 0 (2.1.4) 

Here: u^ and 0^ are the tool surface curvilinear 
coordinates, 4^ is the parameter of motion in the process of 
generation of the gear tooth surface, <|k is the angle of rotation 
of the gear being in mesh with the mating gear. 

Equations (2.1.1) to (2.1.4) provide that the position 
vectors r ^ ^ and r^ 2 ^ and surface unit normals n^ 1 ^ and n^ 2 ^ are 
equal for the gear tooth surfaces in contact (Fig. 2.1). Vector 
equations (2.1.1) and (2.1.2) yield five independent equations 
and the total equation system is 
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fg(u^/6^/i|)^) = 0/ t ^ 2> ^ 2 ^ ~ ^ (2.1.5) 

An instantaneous point contact instead of a line contact is 
guaranteed if the Jacobian differs from zero, i.e. if 

1 2. j 4 5 o / n , 9 i r \ 

D( u ^ , 0 ^ ^ , u 2 , 6 2 , i|j 2 / <j> 2 ) 

If the inequality equation (2.1.6) is observed, then the system 
of equations (2.1.5) may be solved in the neighborhood of the 
contact point by functions 

U l^l^ ,U 2^1^'^l^l^ r /<I > 2^1^ (2.1.7) 

These functions are of class (at least they have continuous 

derivatives of the first order). Functions (2.1.7) and equations 
(2.1.5) enable calculation of the transmission errors (deviation 
of <|> (<|>^) from the prescribed linear function) and the path of 
the contact point over the gear tooth surface. 

For the case when the gear tooth surface is a regular screw 
involute surface, it may be directly represented in a two- 
parametric form and the number of equations in system (2.1.5) may 
be reduced to six. 
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2.2 Simulation of Contact 


Due to the elastic approach of the gear tooth surfaces their 
contact is spread over an elliptical area. It is assumed that 
the magnitude of the elastic approach is known from experiments 
or may be predicted. Knowing in addition the principal 

curvatures and directions for two contacting surfaces at their 
point of contact we may determine the dimensions and orientation 
of the contact ellipse (Litvin, 1968). 

The determination of principal curvature and directions for 
a surface represented in a three-parametric form is a complicated 
computational problem. A substantial simplification of this 
problem may be achieved using the relations between principal 

curvatures and directions, and the parameters of motion for two 

surfaces being in contact at a line. One of the contacting 
surfaces is the tool surface and the other is the generated 

surface . 

Helical gears with modified gear tooth surfaces will be 
designed as surfaces being in point contact at every instant. 
The point of contact traces out on the surface a spatial curve 
(the path of contact) whose location must be controlled. The 


tangent to the path of contact and the derivative of the gear 

ratio (m„, (<{>,)) may be controlled by using the relationship 

d (|) *j^ 2 * J. X 

between principal curvatures and directions for two surfaces that 


are in point contact (ref. 2). Here: 
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is the gear ratio 


2.3 Partial Compensation of Transmission Errors 

Aligned gears transform rotation with a constant gear ratio 
m 2i and 


N 1 

^ 2 0 ^ 1 ^ - N ^ 1 (2.3.1) 

is a linear function. Here: N ^ and H 2 are the numbers of gear 

teeth. An investigation of the effect of helical gear rotational 
axis intersection or crossing indicates that becomes a 

piece-wise function which is nearly linear for each cycle of 
meshing (Fig. 2.2(a)). The transmission errors are determined by 

N 1 

A<f> 2 (<j>^) “ ^2^1^ - ^1 N (2.3.2) 

and they are also represented by a piece-wise linear function 
(Fig. 2.2(b)). Transmission errors of this type cause a 
discontinuity of the gear angular velocity at transfer points and 
vibration becomes inevitable. The new topology of gear tooth 
surfaces proposed in this report allows the absorption of a 
linear function of transmission errors that results in a reduced 
level of vibration. This is based on the possibility to absorb a 
linear function by a parabolic function. 

Consider the interaction of a parabolic function given by 

( 1 ) 2 

A <f> 2 = _a <(> i 


6 

/ 


(2.3.3) 



with a linear function represented by 


A <j> 


( 2 ) 

2 



(2.3.4) 


The resulting function 


A 4> 2 — b<j> - a<j>j 


(2.3.5) 


may be represented in a new coordinate system (Fig. 2.3): 


*2 



(2.3.6) 


where 


2 

*2 = A<t, 2 " 7"2 ' *1 = *1 " 2Z ( 2 . 3 . 7 ) 

4 a 

( 1 ) 2 

Vie consider that a 4> 2 = i- s a predesigned function that 

exists even if misalignments do not appear. The absorption of 

( 2 ) ( 1 ) 2 

function A<t > 2 = b^ by the parabolic function Acf)^ = -a<j>^ means 

that gear misalignment does not change the predesigned parabolic 

function of transmission errors. Thus the resulting function of 

transmission errors A <J> 2 = A^!^ + will keep its shape as a 

parabolic function although the gears are misaligned. The 

resulting function of transmission errors ^ 2 ( ^ 1 ) may be obta i ned 

by translation of the parabolic function a <}> 2 ^ • 

The absorption of a linear function of transmission errors 

by a parabolic function is accompanied by the change of transfer 
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points. The transfer points determine the positions of the gears 

where one pair of teeth is rotating out of mesh and the next pair 

is coming into mesh. The change of transfer points is determined 

h 2 

and A 4>2 = the cycle of meshing of one pair 

2 4a z 

of teeth is <(>. = — i e 1,2. It may happen that the absorption 

1 

of a linear function by a parabolic function is accompanied with 
a change that is too large. If this occurs the transfer points 
and the resulting parabolic function of transmission 
errors, will be represented as a discontinuous function 

for one cycle of meshing (Fig. 2.4). To avoid this, it is 
necessary to limit the tolerances for gear misalignment. 


with Stj = ^ 


2.4 Misalignment of Regular Helical Gears 

Regular helical pinion and gear can be represented by their 
surface position vectors and normal vectors in coordinate system 
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[r 



2 

COS ijl 


n 


sin (<j> p - i|; ) [-t r sing T 


cosi|i„ "■“‘ VM, G '''c' 1 G’ "p 


2 , 

COS l|> 

sin( *G - ^c )[ - t G sinB p 


• 2 

2 Sin 6 d 

t_(cos6 + sin il) 2- 

G p r n cos? 


+ r 2 <j> G J + r 2 sin<j> G 

+ r 2 <(> G ] + r 2 cos <j» G 

r 2^G sin2 ^n tge p 

(2.4.3) 





cos^cos BpCos<}> G + sin^ n sin<f> G 

1! 

CN 

c 

rn„ 1 

2x 

n» 

2y 

= 

-cosi)> n cosBpSin<|> G + simf> n cos<j> G 


n 0 
L 2zJ 


-cos ill sinB 




L n p 


(2.4.4) 


where iji^ and if» c are the pressure ' angles in gear tooth normal 
section and transverse section respectively; B o is the helical 
angle at the pitch cylinder of the pinion and the gear; r^ and r 2 
are the radii pitch cylinder of the pinion and the gear 
respectively; <|>^, and t are the surface parameters of the 
pinion tooth surface; <j> and t_ are the surface parameters of the 
gear tooth surface. 

When the helical pinion and gear are in mesh, with their 
axes misaligned, their position vectors and normal vectors can be 
transformed to the fixed coordinate system S^. The basic ideas 
have already been discussed in the Section 2.1. And the real 
approach of transformation and the matrices to describe the 
misalignment and gear rotation will be given in Section 3.8. 

It is found that when the misalignment occurs the regular 
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pinion and gear surface cannot contact in tangency. That is, in 
Sf, their normals can not be equal in all circumstances. In this 
case, only the gear tooth edge contacts pinion tooth surface in 
tangency. Therefore, in the fixed coordinate system, there are 
four equations to describe the contact, that is, the equalities 
of three position vectors describing gear tooth edge and pinion 
tooth surface as well as the zero product of gear tooth edge 
tangency and pinion tooth normal. 

The computer aided simulation of meshing of misaligned 
helical gears with regular tooth surfaces shows that the 
transformation of rotation is accompanied with large transmission 
errors. There are two sub-cycles of meshing during the complete 
meshing cycle for one pair of teeth. These sub-cycles correspond 
to the meshing of (1) a curve with a surface, and (2) a point 
with the surface. The curve is the involute curve at the edge of 
the tooth of the gear and the point is the tip of the gear tooth 
edge. The transmission errors for the period of a cycle are 
represented by two linear functions (Fig. 2.5). The 
transformation of rotation will be accompanied with a jump of the 
angular velocity of the driven gear and therefore vibrations are 
inevitable . 

The results of computation are presented for the following 
case. Given: number of teeth = 20, N 2 = 40, diametral pitch 
in normal section P n = 10 in - ''', gear tooth length L = 10/P n the 
helical angle 6^ = 15°, the normal pressure angle ^ = 20°. The 
gear axes are crossed and form an angle Ay = 5 arc minutes. The 
computed transmission errors are represented in table 2.4.1. 
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TABLE 2.4.1 - TRANSMISSION ERRORS OF REGULAR HELICAL 
GEARS WITH CROSSED AXES 


*)> 1 / 
deg 

-8 

-5 

CN 

1 

5“ 

T 

7 

mm 

A<f> 

sec 

4.90 

3.06 

1. 22 

-0.61 

-2.45 

-4.29 



2.5 Surface Deviation by the Change of Pinion Lead 


Helical gears in this case are designed as helical gears 
with crossed axes. The crossing angle is chosen with respect to 
the expected tolerances of the gear misalignment (Ay is in the 
range of 10 to 15 arc minutes). The gear ratio for helical gears 
with crossed axes may be represented (Litvin, 1963) as 


M 


12 = 


“2 


r b2 Sin X b2 
r bl Sin X bl 


(2.5.1) 


whe 


re r. i and x, . are the radius of the base cylinder and the 
bi bi 


lead angle on this cylinder, i e 1,2. 


X p2 X pl 


= Ay. 


Here: X . is the lead angle on the pitch cylinder. The advantage 

pi 

of application of crossed helical gears is that the gear ratio is 
not changed by the misalignment (by the change of Ay). The tooth 
surfaces contact each other at a point during meshing. The 
disadvantage of this type of surface deviation is that location 
of the bearing contact of the gears is very sensitive to gear 
misalignment. A slight change of the crossing angle causes 
shifting of the contact to the edge of the tooth (Fig. 2.6). 

The discussed type of surface deviation is reasonable to 
apply for manufacturing of expensive reducers of large dimensions 
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when the lead of the pinion can be adjusted by regrinding. While 
changing by regrinding the parameters r^ and y^ ^ , the 
requirement that the product r b ^ sinX^ must be kept constant. 
Then, the gear ratio M 21 will be of the prescribed value and 
transmission errors caused by the crossing of axes will be zero. 

Theoretically, transmission errors are inevitable if the 
axes of crossed helical gears become intersected. Actually, if 
gear misalignment is of the range of 5 to 10 arc minutes, the 
transmission errors are very small and may be neglected. The 
main problem for this type of misalignment is again the shift of 
the bearing contact to the edge (Fig. 2.6). 

3. GENERATION OF PINION TOOTH SURFACE BY A SURFACE OF REVOLUTION 
3.1 Basic Consideration 

The purpose of this method for deviation of the pinion tooth 
surface is to reduce the sensitivity of the gears to 
misalignment. Also the transmission error must be kept to a low 
level and stabilize the bearing contact. This investigation 
shows that this goal may be achieved by the proposed method of 
crowning but the bearing contact cannot cover the whole 
surface. The reason for this is that the instantaneous contact 
ellipse moves across but not along the surface (Fig. 3.1). 

The proposed method for generation is based on the following 
considerations. It is well known that the generation of a 
helical gear may be performed by an imaginary rack-cutter with 
skew teeth whose normal section represents a regular rack-cutter 
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for spur gears (Fig. 3.2(a)). We may imagine that two generating 
surfaces, E^ and e , are applied to generate the gear tooth 
surface and the pinion tooth surface, respectively (Fig. 
3.2(b)). Surface is a plane (a regular rack-cutter surface), 

VJ 

and E is a cone surface. Surfaces E_ and E are rigidly 
p G p - 

connected and perform translational motion, while the pinion and 
the gear rotate about their axes (Fig. 3.3). The generated 
pinion and gear will be in point contact and transform rotation 
with the prescribed linear function 2 ( ^ l ) * However ' due to gear 
misalignment, function <|> ($ ) becomes a piecewise linear function 
(Fig. 2.2(a)) that is not acceptable. To absorb a linear 

function of transmission errors (Fig. 2.2(b)), a predesigned 
parabolic function of transmission errors is used. For this 
reason a surface of revolution that slightly deviates from the 
cone surface is proposed (Fig. 3.2(c)). The radius of the 
surface of revolution in its axial section determines the level 
of the predesigned parabolic function. The pinion crowning 

process may be accomplished by grinding, shaving or lapping. 

3.2 Principle of Generation and Used Coordinate Systems 

Consider two rigid connecting surfaces E^, and E . The 
generating surface E ^ is a plane and generates the helical gear 
tooth surface that is an involute screw surface. Surface E^ is a 
surface of revolution. Initially, we consider that E is a cone 

hr 

surface and E and E^ contact each other along a straight line 

p G 

that is the generatrix of the cone. 
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Fig. 3.4 shows the generating surfaces e q and E^. Fig. 3.5 
illustrates the process for generation. While the rigidly 
connected generating surfaces perform a translational motion, the 
pinion and gear rotate about their axes 0^ and C> 2 , 
respectively. The parameter of motion of the cutter, S, and the 
angles of rotation of the pinion and the gear, <(> and <J> 2 , are 
related as follows: 


s = r l*p = r 2*G 


(3.2.1) 


where r^ and r 2 are the great centrodes radii, the cutter 

centrode is the straight line that is tangent to the gear 

centrodes. Point I is the instantaneous center of rotation. 

Coordinate systems and S 2 are rigidly connected to the helical 

pinion and the helical gear, whereas coordinate systems S c and S f 

are rigidly connected to the tool surface and fixed frame. The 

generating surface E_ (a plane) and E are covered with a set of 
- G p 

contact lines L G2 and 1^ respectively the instantaneous lines 

of tangency of surfaces E„ and E„ and E and E. (Fig. 2.3, 

Cj / P JL 

a,b). The location of these lines depends on the value of 

parametric <t> G and < f> D ('t > p an< 3 4> G are related). Line L G ^ is the 
line of tangency of generating surfaces E Q and E . 
When E^ and are generated, at any moment, one point on line 
Lp G are generating corresponding point Mp^ and Mq^ on the helical 
gear surface E ^ and E 2 « When the £ and E 2 are meshing without 
misalignment M p j_ and M Gi contact each other in turn. Now it is 

clear why Lp G is not parallel with the edge of E g * The reason is 
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that in this way, the contact ratio of crowned helical pinion and 
regular helical gear will be higher. 


3.3 Tool Surface 


The pinion generating surface is a cone and may be 
represented in an auxiliary coordinate system (Fig. 3.7) as 

follows : 


r 

/V 


d 



u cosQsina 
P 


d-u cosa 
p 


-u sine sina 
P P 


1 


(3.3.1) 


where 0 < u < , 0 < 6 < 2tt . The 

cosa 

represented by 


surface 


normal 


is 


9 £ ( 

N , = ~ 
~d 8 0 


9 £d 

— = v ina 

p 


cosacose 

I 

sina 

cos as in 0 


PJ 


(3.3.2) 


The unit surface normal is (provided u^sina * 0) 


cosacos 0 

P 

n , = sina 
~d 

cosasine 

L pJ 


(3.3.3) 


Figure 3.8 illustrates the installment of the conic tool in 
coordinate system S c step by step as follows: 

(i) From coordinate system to ' , the cone is tangent to the 
plane y^ 1 , o^ ' , ' where the tangent line LpQ (see Section3.2) 
is coincident with y ' axis. Here, we have (Fig. 3.8a). 
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r coso 
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sina 

cosa 

0 
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0 

0 

1 
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-dsina 

dtgasina 

0 

1 


(3.3.4) 


where the d is the height of cone. 


(ii) From coordinate system S ' to S^' the tangent line of cone 
and plane y ^ is declined with an angle y (see Fig. 3.4), also 


the origin o' and o are not coincident, 


Here we have (Fig. 


3.8b). 
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COSlJj 
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(3.3.5) 


where b is the addendum height of the tool. 

(iii) From coordinate system to S a , the tool is declined with 
a pressure angle i|^. For the helical gear, the is measured at 
the normal cross section. Here, we have (Fig. 3.8c). 
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-sinij; 

n 

0 

0 " 




‘"ab 1 

= 

sirnb 

n 

COSlfi 

n 

0 

0 


(3.3. 

6) 



0 

0 


1 

0 






. 0 

0 


0 

1 _ 




(iv) From 

coordinate 

system 

s a 

to 

S c' 

the 

helix angle 

3 P 

considered. 

Here we have (Fig 

. 3. 

8d ) 







"cos 8 

0 

-sin6 

0 " 






p 



P 





[M ] 
ca 

= 

0 

1 

0 


0 


(3.3. 

7) 



sin8 

0 

cos 6 

0 






P 



P 







L 0 

0 

0 


1 . 
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In coordinate system S c , the tool surface and its normal can 
be represented as: 


[r (p) ] = [M ] [M , ] [M, , , ] [M, .] [r , ] 
c ca ab bb 1 bd d 

[n (p) ] = [L ] [L . ] [L, , , ] [L. ,] [n , ] 
c ca ab bb 1 bd d 


(3.3.8) 


where 3 x 3 L matrix is from corresponding 4 x 4 M matrix, 

excluding the last row and last column. Substituting Eg. 

(3.3.1), (3.3.3), (3.3.4), (3.3.5), (3.3.6) and (3.3.7) into 

Eq. (3.3.8), finally we obtain the tool surface E and its normal 

P 

in coordinate system S„ as 


[Pp), 


where 


(P) 

c 

(p) 

c 

(p) 

c 




„'P> 

cx 

n<P> 

cy 

n (p) 

L cz J 


(3.3.9) 


X 


(P) _ 


u cos0 sina (cosacos* cosB + sinasin* cosycosB 
p p y n p y n p 


+ sinasinysinB ) +u sine sina(sin* sinycosB - cosysinB ) 
PPP n M p p 

+ u cosa (cosasin* cosycosB -sinacos* cosB + cosasinysinB ) 
P n p n p p 

d b 

_ ( _ ) [si n* cosycosB + sinysinB ] 

cosa cos^^cosy T n p P 


(P) _ 


u cos0 s ina ( cosas in* - sinacos* cosy) 
p p n y n 


- u sin0 sinacosil siny - ucosa (sinasin* + cosacos* cosy) 
p p y n *n y n 

, , d b 

+ ( - ; ) COS y COS * 

cosa cos* cosy r n 

y n M 
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Z ^ ^ =u cose sina (cosacosib sinB +sinasini|> cosysinB -s inas inycos 8 ) 
c p p r np n p p 

+u sine sina(sinil) sinysinB + cosycosB ) 

P P n p p_ 

+u cosa ( cosas inib cosysinB - sinacosij; sinB - cosa siny) 

Ph k n P n p 


+ ( 


cosa 


COSlb cosy 


) ( -cosysin8_sinib + sinycosB ) 


n 


n 


(P) 


cx 


= cose cosa (cosacosib cosB + sinasinib cosycos B„ 
p n p n p 

+ sinasinysinB ) + sine cosa(sinf sinycosB - cosysinB ) 
p p n p p 

+ sina (sinacosib cosB -cosasinib cosycosB -cosasinysinB ) 
n p n p p 


N 


(P) 

'ey 


= cose cosa ( cosas inib -sinacosib cosy)-sine cosacosiji sinB 
p n n p n p 

+ sina (sinasin^ + cosacos^cosy) 


n^ ^ = cos6pCosa (cos acosi(; n sinBp+s inas in^ n cosys in Bp-sin as inycos B p ) 

+ sine cosa(sim|) sinysinB + cosycosB ) 
p n p p 

+ sina (sinacos^^sinBp -cosas inijj^cosy si nBp+cos as inycos Bp ) 

3.4 Pinion Tooth Surface 

The equation of meshing of the generating surface Ep and the 
helical pinion tooth surface is represented by 


N (P) . V (P1 » = N (P) . (V (P) - V (1) ) = 0 


(3.4.1) 


We can also use equation based on the fact that the contact 
normal of generating and generated surfaces must intersect the 
instantaneous axis of rotation I-I (Litvin, 1968). Thus, we 
obtain (Fig. 2.2). 
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(3.4.2) 


X -x (P) 
c c 


n 


(P) 


cx 


Y -y (P) 
c J c 


n 


(P) 

'ey 


Z 


c 


-z 


(P) 

c 


n 


(P) 

cz 


where (X , Y , Z ) are the coordinates of a point that lies on 

L v 

(p\ (p) 

axis I-I; x , y , z are the coordinates of cone surface; 
c *c c 

n^ v , n„„ and n^„ are the projections of the surface unit 

c x cv v a 

normal. From Fig. 3.5, it is known that 


X = S = rj , Y = 0 
c 1 p c 


Equation (3.4.2), (3.4.3), and (3.3.9) yield 


f,(u ,0 ,d> )=-u [cose ( cosy cos 0 tsinysintii sine ) 
lpp pp P P n P 

+sin9 (sinasinycose -cosacosil) sine -sinasim|j cosysine )] 
P p np n p 

j u 2 2 

+ (—— — — : — r ) [ (cose cos a+sin a). [cospcosB 

cosa cosi^^cosy p p 

+sin^ n sinysinep] -sin6pCosacos^ n sin8p] 

+rd> [cose cosa (cosasintb -sinacosi|> cosy) 

Y p p n n 

+sina (sinasinA +cosacostl) cosy)-sine cosacosij) sintjj ] 

T n n p n n 

= 0 


(3.4.4) 


The helical pinion tooth surface can be obtained by transforming 

the generating tool surface E to coordinate system S-^, together 

P 

with equation of meshing. The coordinate transformation in 
transition from S c to 3 f is represented by the matrix M fc as 
( Fig . 3.5) 

' 1 

[m £ c 1 - o 

0 

. 0 
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1 

0 

0 


0 

0 

1 

0 


-r , <b 
l Y p 

0 

0 

1 


(3.4.5) 


The coordinate transformation in transition from Sj to is 


If 


by the matrix 

Mff as 

(Fig. 

3.5). 



COS * 

p 

sin<t> p 

0 

0 " 



= 

-sin* 

D 

cos* 

p 

0 

0 


( 3.4.6) 


Ir 

0 

tr 

0 

1 

0 




0 

0 

0 

1 . 




Therefore, the helical pinion surface can be represented in 
coordinate system S-^ as: 


“ l M lc llr c P>! = [M l£ 1[M £c llr c P>1 


(3.4.7) 


together with Eq. (3.4.4). Substituting Eq. (3.4.5), (3.4.6) and 

(3.3.9) into Eq. (3.4.7) we have 


[r l ] = 


i p \ ( P ) . 

(x c ir iy cos<t, P +Y c sin4> p 

( p \ ( p ) 

(x c '-ri*p) sin *p+y c cos * p 

*(P) 


(3.4.8) 


Eq. (3.4.8) and (3.4.4) represents the pinion tooth surface 

where y^^ and z^^ are expressed in Eq. (3.3.9). 

c c c 

Using the same approach, the unit normal pinion tooth 

surface can be represented by Eq. (3.4.9) and (3.4.4) 

where n^^, n^^ and n^^ are expressed in Eq. (3.3.9) 
cx cy cz 


[n-jJ 


n cx cos *p + n cy sin *p 

(P) • (P) 

-n sm<f> + n cos* 
cx T p cy p 

(P) 


(3.4.9) 
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It is clear that if we set = O f the pinion becomes spur 
pinion. Therefore, crowning spur pinion using conic tool is only 
a special case of crowning helical pinion. 


3.5 Condition of Pinion Non-Undercutting 

The problem of undercutting of the helical pinion tooth 
surface by crowning is related with the appearance on the pinion 
tooth surface of singular points. From differential geometry, it 
is known that the surface point is singular if the surface normal 
is equal to zero at such a point. 

Litvin, proposed a method to determine a line on the tool 
surface whose points will generate singular points on the surface 
generated by the tool. This line designated by L (Fig. 3.9) must 
be out of the working part of the tool surface to avoid 
undercutting of the pinion by crowning. 

The limiting line L of the tool surface is determined by the 
following equations 


(P) 




( 3. 5. 1) 


f 1 ( V V V = 0 (3.5.2) 


F( v V V = 0 


(3.5.3) 


Vector equation (3.5.1) represents the tool surface [see Eq. 
(3.3.9)]; equation (3.5.2) is the equation of meshing [see Eq. 
(3.4.4)] and equation (3.5.3) comes from the requirement of 
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limiting line L as: 


3 X 


(P) 


3u 


3 Z 


P 

(P) 


3U 


3 f . 


3u 


3 X 


(P) 


39 


3 Z 


p 

(P) 


_c 
'3 9 


3 f 

IF 


V 


(cl) 

ex 


— V 


(cl) 

cz 


3 f 
3 <j> 


1 3 




= 0 


(3.5.4) 


(cl) 

where \r is the relative velocity of the tool and generated 

surface at generating point, represented if coordinate system 

(cl) 

S„. Actually, we can write v as 

C ~c 


v (d) = v (c) _ v (l) 
~c ~c ~c 


(3.5.5) 


where is the velocity of the cutter and is the velocity 

~c ~c 

of the pinion. From Fig. 2.2, we get 


V' 



ds ~ 


_ 


dt 


“p 

= 

0 

= 

0 

i 

0 


. 0 


(3.5.6) 


V ( 1 ) 


~c 


where 


U) 

~P 


U) 

x r (p) +00. 

X UJ 

~P 

~c cl 

~P 

' 0 ' 


r i* P ' 

0 

a) 5 oT = 

-r. 

1 

c 1 

1 



0 


(3. 5.7) 


(3. 5.8) 
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Equations (3.5.6), (3.5.7) and (3.5.8) yield 


v (c!) 


y 


(p) 

c 



+ r, <(> 

1 P 


0 ) 


0 


(3.5.9) 


where X^ and are represented by equations 3.3.9. 

c c 

Equation (3.5.9), (3.3.9), (3.4.3) and (3.5.4) yield 


u YU + 


u P (r i G ' 


+ u * ( WZ + XY) ] + u* XZ = 0 


(3.5.10) 


where 


W 


(B 2 E + A 2 E) + (BD-AF) I + sin 3 e p (A 2 F-B 2 F + 2ABD) 

cos 3 e (2AFB + B 2 D - A 2 D) - sine (2A 2 F + ABD - AEI ) 
P P 

2 

+ cos 9 [BEI + 2B D + ABF] 

P 


2 2 2 
X = [ ( AFs in a - CE)cos0 p + (ADsin a - AEcos a)sin0 p 

+ ( AFcos 2 a - DC)] ( Bcos 0 + Asine + I) 

hr it 

Y = Dtgacos 0 - FtgasinO ^ - Ectgo 

y p P 

Z = cosycosi|> n 
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0* . ( _S 6 ) 

P COSO COSll) cosy 


A = cosycosg + sinili sinysing 
P n p 


B = cosacosil) sing + sinasinil) cosysing - sinasinycosg 
n p r n p p 


C = cosacosib sing 
n p 


D = (cosasiml) - sinacosil) cosy)cosa 
r n r n 


E = (sinasinil) + cosacosil) cosy)sina 
r n v n 


F = cosacosif) siny 


G = Dcose + E - F sine 


I = ctqa (cosasiml) cosysing -sinacosi|) sing -cos as inycos g ) 

n p n p P 


9F 

Using Eq. (3.5.10) and taking it into account that —— * 0, we 

3 6 p 

can represented the as a function of 9 . Undercutting will be 
^ P P 

avoided if 


u (9 ) > 


p p cosa 


(3.5. 11) 


The analysis of Eq. (3.5.10) indicates that the inequality 
(3.5.11) is satisfied if the condition of helical pinion non- 
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undercutting by a regular rack cutter is satisfied. 


3.6 Principal Directions and Curvatures of Tooth Surface 

A simplified approach to determine principal directions and 

curvatures of helical pinion has been proposed by Litvin (Litvin, 

1968). The main idea is representing the principal directions 

and curvatures of generated surface £ ^ by the principal 

directions and curvatures of generating surface £ . 

P 

Let us determine the principal directions of curvatures of 

tool surface £ . The tool surface is a cone surface and its 
P 

principal directions coincide with the direction of the cone 
generatrix and the direction that is perpendicular to the cone 
generatrix . 

The Rodrigue's formula (Eq. 3.6.1) can be used to obtain the 
principal curvature and directions 


K 


I, ii 


v 

~r 


-n 


~r 


(3.6.1) 


where V is the velocity of a point that moves over a surface 
~r 

and n is the derivative of the surface unit normal n, 

~r 

when n changes its direction due to the motion over the surface. 

Using Eq. 3.6.1, the principal directions and curvatures of 
cone surface £ Q can be expressed in coordinate system as: 
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(3.6.2) 


3r, 

3 r . 
~d 


'-sin© ' 
0 

30 

39 


P 

P 


COS0. 


1 


u cos a 
P 


e 


(P) 

II 




cosesina 

-cosa 

sinesina 


1 



0 


The negative sign for indicates that the curvature center is 
located on the negative direction of the surface normal. The 
principal curvatures are invariants with respect to the used 
coordinate system. Whereas the principal direction will be 
represented in coordinate system Sf as 


[e (P) ] 
L §i,n J f 


= t L £c 1 [L cd> 'Sail'd 


(3.6.3) 


where tL cd ] [L ca ] [L ab ] [L ab’ ] [L b’d ] ' 3x3 L matrices can be 
obtained from corresponding 4 x 4 M matrices [see Eq. (3.3.4), 
(3.3.5), (3.3.6) and (3.3.7)], is 3 x 3 unitary matrix (see 
Fig . 3.5). 

The determination of principal curvatures and directions for 
the pinion tooth is based on the following equations (see Litvin, 
1968) 


tan 2 a 


(pl) _ 


2b 13 b 23 


,2 2 (p)_ (pK, 

b 23 b 13 (k I k II b 33 


(3.6.4) 


k (1) -< {1) 


,2 K 2 

b 23 5 13 


_ (ie (p)_ K (p)) b 
' I II '33 


b ^ 3 cos 2cr 


(pl) 


(3.6.5) 
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K (1) + K (1) 

<H -r Kj 


k (p) + Jp) + £l + b \i 

1 11 b 33 


(3.6.6) 


where k^, and e| P ^, e^ are the Principal curvatures and 

unit vector of principal directions of E . and 

p I II 

are the principal curvatures and unit vectors of 
principal directions of E . Angle cr^ P ^ is measured counter 
clockwise from to (Fig. 3.10). The coefficient b 13 , 

b 23 and b 33 have been derived as shown in (Litvin 1968) but 
modified for the case when a rack cutter generates a gear. The 
expressions for b-^ 3 , *^23' an< ^ b 33 are as f°ll° ws: 


b 13 

ii 

r 0 (p) 

e II • £p 


r _ K <p> 

k i 

0 

_ b 23 


-eJP } . u> 

~I ~p 


o 

1 

-< (p) 

ii 


V (pl) e' p) 

v <pl) . e< p) 


(3.6.7) 


b 33 = [n u v[ p) ] + [n ai V (pl) ] - K j p) (V (pl) . ej p) ) 2 
33 ~ ~p ~tr ~ ~p ~ i ~ i 


(p) (pi) (p) 2 

K ll SlI ' 


m 


lp 


m ip 


d4> 


where m. = < m, 

lp ds lp 


dm 


* 


±R 


ds * 


(3.6.8) 


The vectors used in Eq. (3.6.7) and (3.6.8) are represent in 
coordinate system Sf (Fig. 3.5) with i ^ , i^ , k^ as its unit 
vectors of axes. The expressions of the vectors in Eq. (3.6.7) 
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and (3.6.8) are as follows (Fig. 3.5) 


to = a) k _ 
~p ~f 


(3.6.9) 


is the angular velocity of the pinion being in meshing with the 
rack cutter. 


V ( P) 

tr 


-“ r iif 


(3. 6. 10) 


is the transfer velocity of a point on the rack cutter that 
performs translational motion, r^ is the radius of pinion 
centrode . 


PV 


'- (y c +r l> ' 

= a X ri P ^ = at t x f 

= 0J 

x -r. * 
c l y p 

~tr ~p ~f ! f 

L o . 


0 


(3.6.11) 


is the transfer velocity of the pinion 


V (PD = v (p) 



c 

-x +r , d> 
c 1 T 


D 


0 


(3.6.12) 


is the "sliding" velocity - the velocity of a point of rack 
cutter with respect to the same point of pinion. 

Substituting Eq. (3.6.9), (3.6.10), (3.6.11), (3.6.12), 
(3.6.2) and (3.3.9) into Eq. (3.6.7) and (3.6.8), then 
substituting Eq. (3.6.7) and (3.6.8) into (3.6.4), (3.6.5) and 
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(3.6.6), finally we can obtain the principal directions and 

curvatures from Eq. (3.6.4), (3.6.5), (3.6.6) and (3.6.2). Since 
the expressions of tool surface and its principal direction are 
complicated and tedious. It is difficult for us to write down 
the expression of k^, and e^ 2 ^. However, a 

corresponding computer program is developed for determining the 
principal directions and curvatures of any point of the pinion 
surface using the algorithm discussed above. 

The same approach can be used to determine the principal 
directions and curvatures of regular helical gear. However, in 
this case, the generating surface is a plane, and the problem 
becomes simpler. 

3.7 Contact Ellipse and Bearing Contact 

The dimensions and orientation of the instantaneous contact 

ellipse can be determined based on the equations proposed by 

Litvin. (Litvin 1968). The input data for computation is: 

(i = 1,2), a (12) and e , where k! 1 ^ are the principal 

curvature of pinion tooth surface and gear tooth 

surface r is t '^ e u °it vector of the first principal 

( 12 ) 

direction of T. ^ . a is the angle formed by the unit vectors 

of principal directions e|^ and e| 2 ^ (Fig. 3.11) and e is the 
elastic approach of the contacting surfaces. The bearing contact 
is formed by a sot of instantaneous contact ellipses that move 
over the gear tooth surface in the process of meshing. 

The axes of the contact ellipse are directed along 
the r\ axis and e, axis, respectively. The orientation of contact 
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ellipse is determined by the angle a, which is angle 
between n and (Fig* 3.11) The dimensions of ellipse are 
determined by a* and b* which are the half lengths of major and 
minor axis of the ellipse respectively. The following equations 
are used to determine a, a* and b*. 


A = 
B = 


1 

4 

1 

4 


w tan(2aj = 


r (i) (2) |„ 

[K e “ K e " l g l 

r (D ( 2 ) . |„ 

[< £ - < £ + | 9l 


g ^s in2a 
g i~g 2 cos 2cr 


( 12 ) 


( 12 ) 



where 


(3.7.1) 


( 1 ) _ .( 1 ) 


= K ■ 


+ K 


( 2 ) 

II 


( 2 ) _ ( 2 ) 


= < 


+ K 


( 2 ) 

II 


( 1 ) 


- K 


( 1 ) 

II 


g o = k . 


( 2 ) 


- 1C 


( 2 ) 

II 


and 


a* 


e 

i/2 h * _ 

e 

A 

D — 

B 


1/2 


(3.7.2) 


3.8 Simulation of Meshing and Determination of Transmission 

Errors 

The simulation of meshing is a part of the computer aided 
tooth contact analysis ( TCA) program. The simulation of meshing 
is based on equations that provide the continuous tangency of 
contacting surfaces. 

To simulate the meshing of a crowned helical pinion and 
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regular involute helical gear with misaligned axis, we will use 
the following coordinate system, as shown in Fig. 3.11 and Fig. 
3.12: (i) Sj is rigidly connected to the frame (ii) an 
auxiliary coordinate system that is also rigidly connected to 
the frame (iii) and S 2 are rigidly connected to the pinion and 
the gear respectively. The relations between and and 
between S 2 and are shown in Fig. 3.12 and expressed by Mff and 
M h2 as 


M 


f 1 


COS<j> 1 

-sini}) ^ 
0 
0 


s i n <j> 1 

COS<f> ^ 

0 

0 


0 0 ' 

0 0 

1 0 

0 1 . 


(3.8.1) 


M 


h 2 


cosily -sin<()2 0 0 

sin^ cos<(>2 0 0 

0 0 10 

0 0 0 1 . 


(3.8.2) 


It is obvious that the rotation axis of the pinion and gear are 
Z^(Z^) and ^^ 2 ) respectively. Therefore, the error of assembly 
of gears can be simulated with orientation and location of 
coordinate system with respect to Sf. Figure 3.13 shows the 
orientation of S h and S f . When the pinion and gear axes are 
crossed (Fig. 3.13a) and intersected (Fig 3.13b) with operating 
center distance C = O f O h . must be emphasized that C can be 

different from the nominal center distance given by C° = r^ + r 2 , 
where r^ and r 2 are the radii of the pinion and gear pitch 
circles. The coordinate transformation from to Sf is 

represented by [Mf^l as follows (Fig. 3.13) 
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(3.8.3) 


[M 



-cos Ay 
0 

sinAy 

0 


0 sinAy 

-1 0 

0 cosAy 

0 0 


0 ' 
C 
0 
1 


where the gear axes are crossed as shown in Fig. 3.13a. 


t M fhl " 


-1 

0 

0 

0 


0 

-cos A 
-sinAj 
0 


0 

-s inA 
cos 
0 


0 

C 

0 

1 


(3.8.4) 


where the gear axes are intersected as shown in Fig. 3.13b. 

The helical pinion tooth surface and its normal can be 
expressed in coordinate system as: 

[r f 1)] = [M fl ] tr l ] (3.8.5) 

tn f 1)] = [L fl ] [n l ] 

where is 3 x 3 matrix from M^. The helical gear tooth 

surface and its normal can be expressed in coordinate system Sg 
as 

[r< 2) ] = [M £h l[H h2 llr 2 ] (3.8.6) 

* Ihh 1 [L h2 ] 

where and is 3 x 3 matrices from and M h 2 * 

For a helical gear (lefthand), [r2^ anc ^ t n 2^ can 
represented as : 
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[r 2 ] 


2, 

COS <Jj 


cos(^ G -^ n )[-t G sinS p +r 2 ^ G ] + r 2 sin* G 


2 

COS l|) 


^ iT ~- sin(4» G -i| )c )[-t G sine p +r 2 ^ G ] + r 2 cos* G 
I . 2 Sin2f? p . 2 

L t G (cosB p H-sm * n -^)-r 2 * G sin V^p 


[h 2 ] 


— ( c os il> cosg cos(ti^+sin»j) sind>„) 
n p G n G 

cos \b cos& sin<f> -sinij> cosiji n 
n p G n G 


cosili sing 
n p 


(3.8.7) 


where 


and <j> 


are the surface parameters 


and i|> c = arc tg ( tg^ n /cosB p ) is the pressure angle in the 
transverse section of the gear tooth. 

From Eq. (3.3.9), (3.4.4), (3.4.8), (3.8.5), (3.8.6) and 
(3.8.7), it is known that. [substituting Eq. (3.4.4) into(3.3.9) 
to eliminate u p ] : 


r (1) - r (1) U fl a ) 

Ef ” ~f ( VV*1 


ni 1 } = n 1 1 } ( <f> , 9 , <t> i ) 
~r ~r p p l 


( 2 ) _ ( 2 ) . , v 

£f ~ ~f * <t> G f t G , <)) 2 ) 


(3.8.8) 


Ef Ef ^G' t G ,< * > 2^ 


The contact of gear tooth surface is simulated in the TCA program 
by the following equations 
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(3.8.9) 


~f ^p ,0 p ,<j> l^ £f ^G' t G ,<fl 2^ 

£}f ^(♦p»®p» < ^i) = £f ^ = ( $Qt ^g' ^ 2 ^ (3.8.10) 

Vector equation (3.8.9) comes from the equality of the position 
vectors at the contact point and provide three independent 
equations whereas vector equation (3.8.10) comes from the 
equality of surface unit normals and provide only two independent 
equations. Therefore, Eq. (3.8.9) and (3.3.10) yield five 
independent scalar equations as follows: 

f • ( 4> -i , Q = 0 (i = 1,2,. ..,5) (3.8.11) 

11 p p z u o 

Equation system 3.8.11 is the expression in implicit form of 
functions of one variable, i.e. 4 ^. Using theorem of implicit 
function, we can obtain 


*2 = '* , 2 U 1 ) 


(3.8.12) 


considering at the neighborhood of P° = ( <f> 0^, 4)^, 4 * 2 ' t G , 4> G ) where 

P° satisfy Eq.(3.8.11) and 


p (f l,f 2 ,f 3> £ 4 >f 5 ) ^ 

D(<t> p' 9 p' <t> 2 ,t G ,(f> G ) 


(3.8.13) 


The function of transmission error can be determined by 
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(3.8.14) 


N 1 

A<f> 2 ( < f | i ) - “ N - ^1 

3.9 Modification of Generating Surface e^ 

Using cone as a tool surface to generate the helical pinion 
is good way to localize bearing contact. By using the TCA 
program, it can be shown that the transmission errors caused by 
gear misalignment are on a very low level. However, the shape of 
the function of transmission errors is unfavorable as shown in 
Fig. 3.14(a). This shape of the function of transmission errors 
will result in interruption and interference with change of 
meshing gear teeth. A more favorable shape of the function of 
transmission errors is shown in Fig. 3.14(b). To obtain this 
shape of transmission error function, we need to modify 
generating surface E into a revolute surface. The reason for 

p 

using revolute surface is similar to the case for crowning a spur 
pinion (see Litvin and Zhang, 1987). 

The surface of revolution is generated by an arc circle with 
radius p. The arc has a common normal with the cone generatrix 
line at the point M (Fig. 3.15(a)). The circular arc and its 
normal are expressed in an auxiliary coordinate system S 0 as 
follows 


X 

e 


Y 

e 


p[cos(a + 6)-cosct] + ( 


COSot COSll) cosu 
r n 


~ )sina 


“ tLT^' 7 ^T^")sina - 2psin 4 sin(a + 4) 

COSO cost); COSU 2 2 


p [sin ( a+8 ) - sina] + 


cos^cosy 


cos a 


(3.9.1) 
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Id 3 8 

= cosa + 2psin tt cos ( a + ~z) 

cosi|> n cosy 2 2 

Z = 0 
e 


ln e ] 


cos ( a+B ) 
sin( a+B ) 
0 


(3.9.2) 


The surface of revolution is generated by rotation of the 
circular arc about the y e axis and can be represented in 
coordinate system as follows 


(r ] = [L ] [r ] , [n ] = [L ] [r ] 

d de e d de e 


(3.9.3) 


where (Fig. 9.2(b)) 

cos 9 


‘"de 1 ’ 


0 

-sine 


0 sine 

I 

1 0 

0 cose 


Then we obtain 


X d - cosa COSi); cosy 


b )sina-2psin -|sin(a + -|)]cose 


n 


y d = 


h ' ft 6 

+ 2psin tt cos(a + -r) 

cos^cosy 2 2 


Z d = 


(( 


cosa cosil) cosy 
n 


8 B 

)sina-2psin sin(a + -|)]sine 
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(3.9.4) 


‘"a 1 


cos9pCOs(a + 8) 

sin ( a + 8 ) 

sine cos(a + 8) 
P 


The installment of the generating surface in coordinate 
system is the same as described in Section 3.3. The generated 
crowning pinion tooth surface can be obtained using the same 
approach as that described in section 3.4. 

The meshing of gears using the crowning method described in 
this section has been simulated by numerical methods. The 
results of the investigation are illustrated with the following 
example . 

Given: number of pinion teeth = 20, number of gear teeth 
N 2 = 40, diametral pitch in normal section P n = 10 in -1 , pressure 
angle in normal section iji = 20°, helical angle 8 = 15°. The 
pinion tooth is crowned by revolute surface with generatrix 
arc p = 30 in. The revolute surface is deviated from a cone 
(comparing in figure 3.2(b) and (c)). The cone has half apex 
angle a = 20° and bottom radius R = 10.6 in. 

The topology of the pinion tooth surface provides a 
parabolic type of predesigned transmission errors with d = 6 arc 
seconds (Fig. 2.3 (a)) and a path contact that is directed across 
the tooth surface (Fig. 3.1). 

The influence of gear misalignment has been investigated 
with the developed computer program and the results of 
computation are represented in table 3.9.1 and 3.9.2 for crossed 
and intersected gear axes, respectively. The misalignment of 
gear axes is 5 arc minutes. 
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The results of computation show that the resulting function 
of transmission errors is a parabolic one. Thus the linear 
function of transmission errors that was caused by gear 
misalignment has been absorbed by the predesigned parabolic 
function. Fig. 3.14 show the results of transmissions errors for 
crossed helical gears with (Table 3.9.1) and without (Table 
3.4.1) crowning of the pinion. 

TABLE 3.9.1 - TRANSMISSION ERRORS OF CROSSED HELICAL GEARS 


(deg) 

-14 

-11 

— 

oo 

' 

-5 

-2 

1 

4 

A<|> ^ 
(sec ) 

-4.99 

-1.51 

0.65 

1.51 

1.05 

-0.75 

-3.87 


TABLE 3.9.2 - TRANSMISSION ERRORS OF INTERSECTED HELICAL GEARS 


tieg) 

-11 

i 

CO 



-5 

-2 

i— 1 

4 

7 

(sec ) 

g 

-2.72 

-0.60 

0.20 

-0.32 


-5.40 


4. CROWNED HELICAL PINION WITH LONGITUDINAL PATH CONTACT 
4.1 Basic Concepts and Considerations 

A longitudinal path of contact means that the gear tooth 
surfaces are in contact at a point at every instant and the 
instantaneous contact ellipse moves along but not across the 
surface (Fig. 4.1(a)). It can be expected that this type of 

contact provides improved conditions of lubrication. Until now 
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only the Novikov-Wildhaber ' s gears could provide a longitudinal 
path of contact. A disadvantage of this type of gearing is their 
sensitivity to the change of the center distance and the axes 
misalignment. The sensitivity to non-ideal orientation of the 
meshing gears cause a higher level of gear noise in comparison 
with regular involute helical gears. Litvin et al. (Litvin, 
1985) proposed a compromising type of non-conf ormal helical gears 
that may be placed between regular helical gears and Novikov- 
Wildhaber helical gears. The gears of the proposed gear train 
are the combination of a regular involute helical gear and a 
specially crowned helical pinion. The investigation of 
transmission errors for helical gears with a longitudinal path of 
contact shows that their good bearing contact is accompanied with 
an undesirable increased level of linear transmission errors. 
The authors propose to compensate this disadvantage by a 
predesigned parabolic function of transmission errors, that will 
absorb the linear function of transmission errors (see section 
2.3). The two following methods for derivation of the pinion 
tooth surface will the modified topology will now be considered. 


4.2 Method 1 

Consider that two rigidly connected generating 

surfaces, and e , are used for the generation of the gear and 
G p 

the pinion, respectively (Fig. 4.1(b)). Surface e q is a plane 

and represents the surface of a regular rack-cutter; 

surface E is a cylindrical surface whose cross-section is a 
P 

circular arc. We may imagine that while surfaces E and E 
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translate, as the pinion and the gear rotate about their axes. 
To provide the predesigned parabolic function of transmission 
errors it is necessary to observe the following transmission 
functions by generation 

N 

— = r „ = const, — = r 0 (~ - 2a<$> ) = f ( <J> ) (4.2.1) 

o) G 2 “p 2 n 2 P P 

Here: u> G and are the angular velocities of pinion and gear by 

cutting; V is the velocity of the rack-cutter in translational 
motion; N-i and No are the gear and pinion tooth numbers; <)> is 

j. z p 

the angle of rotation of the pinion by cutting. .The generated 
gears will be in point contact at every instant and transform 
rotation with the function 

N 

4*2^1^ = ^1 ** a< ^l 0 < <() ^ < jp- (4.2.2) 

This function relates the angles of rotation of the pinion and 
the gear, ^ and $ respectively, for one cycle of meshing. 

The predesigned function of transmission errors is 


**2 



(4.2.3) 


It is evident that after differentiation of function (4.2.2) we 
obtain that the gear ratio u> 2 /u)^ satisfies equation (4.2.1), 
if 4 > and <|> 2 are used instead of and <j> G . 

To apply this method of generation in practice it is 
necessary to vary the angular velocity of the pinion in the 
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process of its generation that may be accomplished by a computer 

controlled machine for cutting. 

It is obvious that in this method the gear tooth surface is 

kept as regular skew involute surface as shown in Eg. (2. 4.1) and 

(2.4.2) since its generating surface £ is a plane and generating 

G 

motion is V = go r . Now, let us derive the pinion tooth 

G 2 

surface equations. 

The pinion generating surface is a cylindrical surface 
with circular arc as its cross section and can be represented in 
an auxiliary coordinate system S as follows: (see Fig. 4.2) 


r = y a 
~a 

z 

a 


R [cos ( <j> n + a ) - cos^] 
R[sin(ip n +a) - sin^] 


(4.2.4) 


The unit surface normal is represented by 


n = n 
~a ay 


cos ( ^ n +oc ) 

= sin(^i +a) 
n 


(4.2.5) 


From coordinate system S a to S c the helical angle 8 is 
introduced. Based on Fig. 3.8(d) and [M ca ] as shown in Eq. 
(3.3.7), we have 


[M ] r 
ca ~a 


R [cos (<!> n +a) 
R[sin(iJ) +a ) 
R [cos ( 'l' n + a ) 


- cost)) ] - t sin 8 

n p p 

- sin^ ] 

- cosifi ] + t sing 

11 ir 


(4.2.6) 
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n 


(P) 

cx 


n 


(P) 

cy 

(P) 


cz 


[L ]r 
ca ~a 


cos ( f +a)cosB 
n p 

sin(<{) +a ) 

cos(ii> +a)sing 
n p 


(4.2.7) 


The generating process is described in the Fig. 3.5 where 


N. 


s ' r 2*p ( ii; - a V 


(4.2.3) 


The equation of meshing of the generating surface z^ and the 
helical pinion tooth surface is given by: 


(P) 


V 


(PD _ Jp> 


= n 


(V 


(P) 


- v 11 ’) 


(4.2.9) 


In the system shown in Fig. 3.5, it can be found that 
ds 


V 


(P) _ 


dt 

0 

0 J 


V 


( 1 ) 


-r , -y 
1 1 c 


-S + X 


(P) 

(P) 


(4.2.10) 


d (j) 

Taking into account that -rr-2- = w and from Eq. (4.2.1) we have 

dtp 


V v 


(pi) = v (p) _ v (l) = 


V' 


-r 2 (jji - 2ay + rj+y^ 1 


I \ N 1 

-X ^ +r 0 <j> (==■ - a 4 > ) 
c 2 y p N 2 *p 


0 
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+ 2r~ad> 
2 


(P) 


, \ N, 

-x p + r 9 (f > (— 
c 2 Y p'N 2 

0 


ad> ) 
P 


(4.2.11) 


Substituting Eq. (4.2.6), (4.2.7) and (4.2.11) into Eq. (4.2.9). 

the equation of meshing can be obtained as: 


cos(<|) +a)cos8 { 
n p 1 


R[sin(iJ) +a)-sintj) ] 
n n 


+ 2r a<J> } 

2 p - 


+ sin(<J) n +a){ -R [cos ( i|; n +a ) -cost); n ] 
= 0 


N 1 

tsin 6p | + r 2 *p ( N; ' a V 

(4.2.12) 


Equation of meshing gives the relation among three variables, 

that is, t , a and <|> . It is obvious that from Eq. 4.2.12, the 
P p 

equation of meshing can be rewritten as: 


t = cotan(^ +a)cotan(/? ) { R [sin(i|> +a)-sin 4 » ]+2r 9 a<j> } 
p ii p n 1 1 « y 


i N i 

+ iT^T {-R[cos(^ n +a)-cos4- n ] + " a V } 


(4.2.13) 


The pinion tooth surface can be determined with the same approach 
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as that described in the section 3.4 and represented by 
Eq. (3.4.8) and (3.4.9) together with Eq . ( 4 . 2 . 6 ) , ( 4 . 2 . 7 ) and 
(4.2.13). 

4.3 Method 2 

The derivation of the crowned pinion tooth surface is based 
on two stages of synthesis. On the first stage it is assumed the 
pinion tooth surface z^ is exact conjugate surface to the gear 
tooth surface which is regular skew involute surface under the 
condition that the rotation transformed by the pinion and the 
year is described in Eq. (4.2.2) with predesigned parabolic 
function of transmission errors. 

On the second stage of synthesis it is necessary to localize 
the bearing contact and substitute the instantaneous line contact 
by a point contact. This becomes possible if the pinion tooth 
surface will be deviated as it is shown in Fig. 4.1(c). Only a 
narrow strip, L, will be kept while z ^ will be changed into zj. 
The deviation of z^ with respect to Z ^ may be accomplished in 
various ways, for instance, in such a way, that the cross-section 
of z^ is only a circular arc. The generation of zj requires a 
computer controlled machine to relate the motions of the tool 
surface and the being generated pinion surface z|. The tool 
surface (it may be only a plane) and z^ will be in point contact 
in the process of generation. 

Now, let us derive the pinion tooth surface equations based 
on the two stages discussed above. First, we can derive the 
intermediate pinion tooth surface z^ as that generated by gear 
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tooth surface E ^ . The relation of rotation for the pinion and 
the gear is shown in Eq.(4.2.2). The gear tooth surface e 2 and 
its unit normal are shown in Eq. (2.4.3) and (2.4.4). 

To represent the gear tooth surface e 2 in the fixed 
coordinate system Sf in Fig. 4.3, it is necessary to transform 
[r 2 ] and [ n 2 ] into [r^ ] and [n^ ] as: 



[M f2 l [ r 2 ] 





tn< 2) l • 

[L f 2 J t n 2 l 




(4.3.1) 


-COS ({) 2 

sine)) 2 

0 

o' 


[M £2 ) - 

-sirup 

-cos<)> 2 

0 

0 



0 

0 

1 

c 



0 

0 

0 

1 . 



Substituting Eq. (2.4.3) and(2.4.4) into Eq. (4.3.1) we obtain: 


[r< 2) ] = 


,( 2 ) 

‘f 

I2! 


2 , 

COS to 


cos^ c cos(t G ^ c -^ 2 )[-t G sin 6 p +r 2 4 , G ]-r 2 sin(^ G - 4 . 2 ) 

2 

Cos 1 1> 

C - C oTT" Cos( ^G^C -♦ 2 )[-t G sinB +r 2 * G ]-r 2 sin(* G -* 2 ) 


2 2 2 

t G cos6 p ( 1+sin i|i n +g 8 p )-r 2 <(> G sin <(/ n tgB p 


(4.3.2) 



r n (2) 1 

n fx 


cosiJ) n cos 8 p cos ( <t> G “4> 2 ) +simj» n sin ( <f> G “'(> 2 ) " 

fn (2) l - 

[n f J = 

n (2) 

fy 

= 

-cosij) n cos8 p sin(((i G -<t> 2 )+sinij) n cos ( 4> G — 4> 2 ) 


[ n fz J 


cosii) sins 

L n p J 
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The generating process is shown in Fig. 4.3 
where 4> and 4> ^ is given by Eg. (4.2.2). The equation of meshing 
of the generating surface T. ^ and generated surface E ^ is 
described as: 


n ( 2) . v (21 > . N (2) (V (2) - v (1) ) - 0 


(4.3.4) 


( 21 ) 

The relative velocity V can be expressed as 


y( 21 ) = V (2) - V (1) = u>~ x '(0,0. + r f ) - x r f 

/V /V /V ^ ^ | /v /v, 3l 


where 


(4.3.5) 


( 1 ) . 


0 

0 

•1 


a) . ; to 
1 ~ 


( 2 ) _ 


0 

0 

1 


0) n = 


0 

0 

1 


N 1 

“l [ N 2 “ 2a + l J 


° 2 ° 1 = 


0 

-C 

0 


( 2 ) 

f 

( 2 ) 

f 

( 2 ) 

f 


and the relation between and to 2 is obtained by taking 
derivative of Eg. (4.2.2). By rearranging and simplifying Eq. 

(4.3.5), we obtain 


N- 



r v (2)i 




- Y f 

N 2 

o)-^ [jj— + 1 — 2a 2 ^ + 

r c i 

V (21) = 

CM 

1 ^ U-l 

X o 

i 

» 

o o 


’i [ n 2 “ 2a<t> l ] 


(4.3.6) 


Substituting Eg. (4.3.6) (4.3.2) and (4.3.3) into Eq. (4.3.4) and 
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simplifying it, we can obtain: 


2a<j> . 


cos ( <(> g ~4> 2 ~^ c ) cosi|» c [l n 1 /N 2 


(4.3.7) 


where i| » is the pressure angle in transverse section of helical 
gear. Equation 4.3.7 is the equation of meshing which describes 
the relation between <j> 0 and <|> 0 . We can rearrange Eq. 4.3.7 as: 


N 1 2 

= ^2 + *^1^ = N - ^l -a ^l + ^ ^ 4* i ^ 


(4.3.8) 


2a 


where X(((> 1 ) = i|> c -arccos cos^Ll- — ] 

The pinion tooth surface i ^ and its normal can be determined 


( 2 ) 


by transforming r^ into coordinate system S-^ as 


[r x ] = [M 

If 1 [ E£ 2>1 




(4.3.9) 

i-j 

ii 

c 






COS((> 1 

-sin<|> 

0 

0 ' 

< 

where [M.^] = 

sin<j> ^ 

cos<f> ^ 

0 

0 


0 

0 

1 

0 



0 

0 

0 

1 



Substituting (4.3.2) (4.3.3) into Eq. (4.3.9), we can get 
x. 


[r,] - 


[n x ] = 


n 


lx 

l iy 


n 


lz 


(4.3.10) 
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where 


2 , 

COS lp 


n 


X 1 = cos(4» 1 +4' c -X(<|) 1 ))[-t G sinB p +r 2 ^ G ] 

+ r 2 sin { 4> ( <p ) ) - Csin^ 

2 

COS lp 

y l = ^sT c n ' sin (^i + 'l'c" X( h )[ " t G sin6 p +r 2^G ] 

- r 2 cos ( <p -X ( <p ^ ) ) + Ccosip^ 


2 2 2 

Z 1 = t Q cosf? (1+sin (P n +g 3 ) - r 2 <p G sin ip^t 6 


g P 


lx 


cosip cos$ 

_n £ 

COSlp 


cos ( <p 1 — X ( <p x ) +ip ) 


COSlp cosS 

n ly = 2 sinUj-xUjMi.p 


n. = cosip sin3 
lz n p 


Equations (4.3.10), (4.3.8) and (4.2.2) represent pinion 
tooth surface E^ in coordinate system S^. But E ^ is only an 
intermediate surface. Now we come to the second stage, that is, 


to deviate E-^ into E j . 

The 

surface e^ must 

satisfy 

such 

condition that when 

it is 

in 

mesh 

with 

gear 

tooth 

surface E 2 without misalignment, 

the 

contact 

path will 

be in 

longitudinal direction. 

The surface 

1 . , , 
E j will 

be 

formed 

in two 


steps: (i) the contact path curve in E ^ is chosen and kept, 

(ii) the profile of surface E ^ in transverse section is replaced 


48 


by a smaller circular arc attached on the L in such a way that 
the original normal of E^ along L is kept. 

To choose curve L on E ^ we could define that the contact 
path on the gear tooth surface is in the middle of the tooth, 
that is 


r 



(4.3.11) 


where r = /x 2 + y 2 is the distance of the surface point to gear 
axis in transverse section. Substituting Eq. (2.4.3) into Eq. 
(4.3.10) and simplyfing it, we can obtain the relation of surface 
parameter along the contact path as 


C 2+G ' t G slnS p = ° 

Applying Eq. (4.3.12) to Eq. 
pinion tooth surface as: 

x 1 = r 2 sin( 4 > 1 -A( 4 > 1 ) ) - 

y l = _r 2 cos i -x i )) 
N i 

Z 1 = r 2 ct 9 B p ( N; *1 ' 3 
where the <j> 1 is the curve 


(4.3.12) 

(4.3.10), we can find the curve L on 

Csin<(> i 
+ Ccos<f)^ 

$1 + A(* 1 )) (4.3.13) 

parameter and X ( c(> ^ ) is described in 


Eg. (4.3.8). 
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Now we should attach the circular arc to the L described by 
Eq. (4.3.13) in the transverse section to form new surface e*. 
Also the tangent of the arc at the point on the line L must be 
perpendicular to the normal of £ ^ described in Eq. (4.3.10). As 
shown in Fig. 4.4, the equation of new surface E ^ are: 

x^ = x^ + r [cos ( y +a ) -cosy ] 

y^ = y^ + R [sin ( y +a ) -siny ] (4.3.14) 


where y = <j> , - A ( $ . ) +ip . Equation (4.3.14) describes the new 
X 1 c 

pinion surface E^. In Eq. (4.3.14) when a = 0, the designed 
contact path L is obtained. 


4.4 Discussion and Example 

In section 4.2 and 4.3, two methods have been presented with 
derivation of the equations of pinion tooth surface. Comparing 
the two methods for the generation of the pinion tooth surface, 
it may be concluded that both provide a localized bearing 
contact, a longitudinal path of contact and predesigned parabolic 
function of transmission errors. The difference between these 
methods is that the tool and pinion tooth surfaces are in line 
contact by applying the first method for generation and in point 
contact by the second one. The disadvantage of both methods for 
crowning of the pinion is that the transmission errors caused by 
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gear misalignment are large and it is necessary to envision a 
high level of the predesigned parabolic function for the 
absorption of transmission errors. This is illustrated with the 
following example (the algorithms for simulation have been 
discussed in Section 3.8). 

Given (the data is from ref. 3): pinion tooth number N-^ = 

12, gear tooth number N 2 = 94; diametral pitch in normal section 
P n = 2 in - -*-; pressure angle in normal section iji = 30°; helical 
angle 0 = 15° . 

The pinion tooth surface is a crowned surface whose cross- 
section is an arc of a circle of the radius 0.3584in. The 
predesigned parabolic function is of the level d = 25 arc seconds 
(Fig. 2.3(a)). 

Consider now that the axes of the gear and the pinion are 
crossed and the crossing angle is 3 arc minutes. The computer 
program for the simulation of meshing provides the data of 

transmission errors that is given in Table 4.1. The data of 

Table 4.1 shows that the resulting function of transmission 
errors is a parabolic function. Thus, the linear function of 

transmission errors caused by misalignment of gear axes has been 
absorbed by the predesigned parabolic function. 

Table 4.2 represents the transmission errors for the same 
helical gears for the case when the gear axes are intersected and 
form an angle of 3 arc minutes. The resulting function of 

transmission errors is again a parabolic function with the level 
d = 26.2 arc seconds. The relatively high level of transmission 
errors is the price that must be paid for the longitudinal path 
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of contact. However, the proposed topology of the pinion tooth 
surface provides a reduction of the level of gear noise since the 
linear function of transmission errors is substituted by a 
parabolic function. 

TABLE 4.1 TRANSMISSION ERRORS FOR CROSSED HELICAL GEARS 




<f> 2 

(deg ) 

-23 

-18 

-13 

-8 

-3 

2 

7 

(sec ) 

-17.94 

-3.06 

5.64 

8.23 

4.84 

-4.39 

-19.37 


TABLE 4.2 TRANSMISSION ERRORS OF INTERSECTED HELICAL GEARS 


^ 1 

(deg) 

-20 

-15 

-10 

-5 

0 

5 

10 


-23.06 

-8.50 

0.15 

2.96 

0.00 

-8.66 

-22.95 


5. DEFORMATION OF HELICAL GEAR SHAFT 
5.1 Basic Concepts and Considerations 

Deformation of gear shafts always exists when the gears are 
used to transmit power. This is because under the load the 
force applied on gear tooth surface is transferred to the gear 
shaft and the shaft is not rigid body. It can be proven that the 
deformation of gear shaft results in the same effects as 
misalignment induced by assembly. The misalignment from assembly 
could be reduced to as little as possible. But the shaft 
deformation is inevitable. Fortunately, all the transmission 


52 



















errors and shift of bearing contact due to deformation of gear 
shaft can be compensated by crowning helical pinion tooth surface 
with the methods discussed in Chapters 3 and 4. 


5.2 Force Applied on Gear Shaft 

When the pinion and gear mesh a force, between their tooth 
surfaces, is applied on the contact point. The direction of the 
force is along the normal of the contact point. For the case of 
regular helical gears, in the fixed coordinate system Sj as 
described in section 3.8, the force direction is a constant since 
the normal of the contact point is a constant (Litvin, 1968) and 
can be expressed as : 


n = 
~c 


cosib cosg 1 
r n p 1 


sim|> 


n 


cosip sine 
n p 


(5.2.1) 


As shown in Fig. 5.1, the force applied on the pinion tooth 
surface is: 


F = F, + F = -F 
~ ~t ~z 


COS>p 

n 

simli 

n 

COSll) 

n 


cos S I 

P 


sing 

P J 


(5.2.2) 



cosib cosB^ 
n p 


0 

where F^ = -F 

sirnl; 

n 

; F = -F 
z 

0 


0 

L 


cosi|) sing 

L n PJ 
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F is called as transverse force since it is applied on the 

transverse section of the pinion and the gear and F^ is called as 

axial force. Transverse force F, is always in tangency with the 

~t 

base circle in tranverse section. Therefore, if the torque 
transferred by the gears is constant, the magnitude of transverse 
force is also constant. So is the magnitude of axial force since 
it has certain ratio with transverse force. 

Both transverse force and axial force can be transferred to 
the axis of gear shaft with resultant 

torque (designated by F* and F*) . For the transverse force, the 
resultant torque is the torque transferred by the gear. For the 
axial force, the resultant torque is balanced by the support 
bearings. Assuming the force is applied on the middle section of 
the gear, after the force is decomposed and transferred to the 
axis of gear shaft as shown in Fig. 5.1 (b), both the transverse 
force and axial force will act on the origin of coordinate system 
Sf. Actually, both forces will cause deformation of the shaft. 
But since the axial force only cause very small tension or 
compression of the shaft, it can be neglected. 

5.3 Modelling of Shaft Deformation 

As shown in Fig. 5.2, the transverse force F is applied on 
the point A which is the center of the shaft corresponding to the 
pinion or gear middle cross section. Under the force, the 

deformation of the shaft is composed of two parts, that is, 
deflection of shaft at the point A designated by Vp and rotation 
of the shaft cross section with A as a center designated by 


54 



A . The value and A depend on the magnitude of transverse 
P P P 

force, geometry and material of shaft and the way how the shaft 
is supported. (Timoshenko 1973). 

Now, let us model the deformation of the helical pinion 
shaft. For convenience and consistency with the previous 
chapters, We establish our coordinate system as follows: 

(1) As shown in Fig. 5.3a 3^ is a fixed coordinate system and 

rigidly connected to the frame. S a is an auxiliary coordinate 

system with the Y q axis coincident with the direction of 

transverse force. The angle i|j is the pressure angle in the 

transverse section of helical gear. The matrix [Mf a ] can 

transfer vector from S Q to Sf and is expressed as: 

sinii; cosip 0 0 

c T c 

-cos>|) simfi 0 0 

c c 

o o 10 

o o 0 1 

(2) Coordinate system S fl and S Q , [Fig. 5.3(b)] are connected to 
the axis of pinion shaft. Without deformation, Sf ' and S a ' are 
coincident with Sf and S a respectively. The matrix [M a * f 1 3 can 
transfer vector from S f ' to S a i and is expressed as: 

simp -cosij> 0 0 

c c 

cosil> sin\J> 0 0 (5.3.2) 

c c 

0 0 10 

0 0 0 1 _ 

(3) Coordinate systems S a and S a * are not coincident when shaft 
deformation occurs. As shown in Fig. 5.2, their relation can be 
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expressed as: 


M 


aa 1 


1 

0 

0 

0 


0 

cos X 

I 

sinX 


0 

-s inX 

I 

cos X 


0 I 

“V ! 

P 

0 I 

1 j 


(5.3.3) 


(4) As described in previous chapter, pinion tooth surface can 
be expressed in coordinate systems S ^ as a position 
vector r 1 with its normal n . The relation of S-,^ and S f ' is 
shown in Fig. 5.3(c) and expressed as: 


M 


f ' 1 


COS<|> 1 


0 

0 


sin*. 0 0 

cos^ 0 0 

0 1 0 

0 0 1 


(5.3.4) 


After the coordinate systems are established, it is easy to find 
that the deformation of the pinion shaft can be modelled by 
matrix rtif f ’ as: 


Wff.J = tM f J W.,,1 


fa 


aa 


a'f 


1+cos * (cosX -1) ,sin* cos* (cosX -1) , -cos* sinX 
c p c c p c p 

2 

sin* cos* (cosX -1) , 1+sin * (cosX -1) , -sin* sinX 

'-up C p C p 


sinX cos* 
p r c 


sinX sin* 
p v c 


cos X 


, — V cos * 
p c 

, -v sin* 
p r c 

, 0 
1 


(5.3.5) 
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Since X is very small angle, it is reasonable to use X instead 
p P 

of sinX and one instead of cosX . Therefore, Eq. (5.3.5) can be 
P P 

written as: 


[H££,] 


l 

0 

X cosi|> , 
P c 

0 


0 

1 

X sintl) , 
p c 

0 


■X COSlb 
P C 

■X sinij* , 
P c 

1 

0 


-V cos ili 
p c 

-V s i n il> 

P c 

0 


(5.3.6) 


Also, for the gear shaft, the coordinate systems S Q , S G ,, 




and 


V 


are used instead of 


>f ' 


S f \ 


S. 


and 


V' 


and X^, and V_ are used instead of X and V . Then, using the 
G G P P 

same ideas, the deformation of gear shaft can be modelled by 
matrix [m ' ] as: 


M 


GG' 


1+cos i|> c (cosX G ~l ) , sim|) c cost() c (cosX G ~l) , cosi|> c sinX G , V G cos^ c 

2 

sini|> cosij) (cosX G -l), 1+sin ij> c ( cos X G ~1 ) , cos^sinXg , V G cos^ c 


-sinX^cosib 
G T c 


, -sinX^sindi 
G c 


cos X 
0 


0 

1 


Or considering that X„ is very small angle: 


(5.3.7) 


[M GG ' 1 


X^COSll) 

G r c 


0 

1 

■X_sint|) 

G T c 


+ X~COSiJj 

G c 

+ X„sint|> 

G c 

1 

0 


V G cos *c 

V G sin *c 

0 

1 


(5.3.8) 


It must be emphasized that the gear tooth surfaces are expressed 
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in coordinate system S 2 . The relation of S 2 and Sg 1 
Fig. 5.4 and expressed as: 


is shown in 


[M 


G' 2 J 


-cos <)>2 sin<|> 2 0 0 

-sin<f>2 -cos<|>2 0 0 

0 0 10 

0 0 0 1 


(5.3.9) 


Also as shown in Fig. 5.5, the coordinate system Sg and Sq are 
not coincident. There is a center distance C between their 
origins. Therefore is represented as: 


[M 


f G 


10 0 0 
0 1 0 C 

0 0 10 

0 0 0 1 


(5.3.10) 


Now it is easy to simulate the performance of the gears with the 

deformation of their shafts. Assuming, there is a helical pinion 

in coordinate system represented by position vector r^ and 

normal vector n^ and a helical gear in coordinate system S 2 

represented by r„ and n 0 , we can write the tooth contact equation 
~ Z ~ Z 

as : 


Wff 

J 

‘W 

'El 1 ' 

‘"fG 1 

[m gg ' ] 

[M G' 2 ] 

'E 2 1 

[L £f 

,] 

IL f -l 1 

'El 1 ' 

IL fG ] 

[ l gg ' ] 

[L G' 2 ] 

t S.2 ] 


where L matrices are 3x3 matrices from corresponding 4 x 4 M 
matrices, crossing 4th row and 4th column and M matrices can be 
found in this section. Applying the same ideas discussed in 
Chapter 2 and used in section 3.8, we can find the transmission 
errors and the shift of the bearing contact by computer aided 
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simulation. 


It is interesting to mention that the deformation of the 
shaft can be expressed as a combination of misalignment with 
crossing axes, misalignment with intersecting axes, and change of 
the center distance. This is because, for example 



1 

0 

-X COSlJi 

p c 

-V cos ib 
P 

tM f f ' J 

0 

1 

-X sinib 
p *c 

-V s i n ib 
P 


X COSlii 

p r c 

X sinib 
p r c 

1 

0 


0 

0 

0 

1 



" 

1 


0 

-X COS lb 0 


' 1 0 

0 

0 ~ 






P 

c 





= 


0 


1 

0 

0 


0 1 

-X sinib 

0 










P c 



X 

COSlb 

0 

1 

0 


0 X sinib 

1 

0 


P 


' c 





p v c 




- 

0 


0 

0 

1 


_ 0 0 

0 

1 

1 

0 


0 


-V cosiJj 










p r c 





0 

1 


0 


-V si nip 



(5.3. 

12) 






P c 





0 

0 


1 


0 





0 

0 


0 


1 





where the 

three 

decomposed matrices 

are represented 

by the matrix 

for 

crossing 

axis with 

small 

angle X cosip 
P c 

, matrix 

for 


intersecting axis with small angle X sin^, and matrix for axis 


displacement . 


5.4 Example and Discussion 

The results of investigation in this chapter are illustrated 
with the following example. Given: number of pinion tooth 
N^= 20, number of gear tooth N 2 = 40, diametral pitch in normal 
section P n = 10 in -1 , pressure angle in normal 


59 


section = 20° , helical angle 3 = 15°, gear tooth length L = 
5/P n . Also assume deformation values Vp = Vq = 0.0125 in, 

A = A = 2 min. The computer program for simulation provides 

ir ' J 

the data of transmission errors in Table 5.1 

From Table 5.1, it is known that the transmission errors due 
to gear shaft deformation is an approximately linear function for 
regular skew involute pinion and gear. Similar to the case of 
gear axis misalignment, the linear function of transmission 
errors can be absorbed by predesigned parabolic function of 
transmission errors that is obtained by crowning helical pinion 
tooth surface. 

TABLE 5.1 TRANSMISSION ERRORS FOR HELICAL GEAR WITH 
DEFORMED SHAFT 


CT 

Q) 

r^O 

3 

6 

9 

12 

15 

8 

21 

a<|> 9 
( sec ) 

1.38 

3. 16 

4.74 

6.32 

7.90 

9.48 

11.06 


6. CONCLUSION 

Several methods of crowning helical pinion tooth surface 
have been developed. The modified pinion tooth surface can 
provide predesigned parabolic function of transmission errors 
that are able to absorb linear function of transmission errors 
induced by misalignment. Also, the modified pinion tooth surface 
can improve the bearing contact. Principles of computer aided 
simulation of meshing, contact, and respective computer programs 
have also been developed. The numerical results of examples of 
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cronwed helical pinion in mesh with regular helical gear show 
that the ideas of crowning are useful to get favourable bearing 
contact and allowable transmission errors. But the synthesis of 
pinion tooth surface should be based on a compromise between the 
requirements of transmission errors and the patterns of the 
bearing contact. 
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FIG. 3.2 Generating Surfaces 
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FIG. 3.3 Generation of Gears 
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FIG. 3.5 Coordinate Systems for Generating Pinion and 
Gear Simultaneously 
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FIG. 3.6 Pinion and Gear Generating Surfaces 
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FIG. 3.10 Principal Directions of Generated and Generating Surfaces 
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FIG. 4.3 Basic Coordinate Systems for Meshing Gears 
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FIG. 5.2 Shaft Deflection 
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o o o n n 


******************************************************* 


C. . 
C. . 

c. . 
c. . 
c. . 
c. . 
c. . 
c. . 


* * 

* PROGRAM I * 

* SURFACE OF HELICAL PINION GENERATED BY CONE CUTTER * 

* * 

* AUTHORS: FAYDOR LITVIN * 

* JIAO ZHANG * 

* * 


C. . . * * 

£ ******************************************************* 

c 

C PURPOSE 
C 

C THIS PROGRAM IS USED TO CALCULATE THE SURFACE OF A HELICAL PINION 
C WHICH IS GENERATED BY CONE CUTTER 
C 

C NOTE 
C 

C THIS PROGRAME IS WRITTEN IN FORTRAN 77. IT CAN BE COMPILED BY V 
C COMPILE IN IBM MAINFRAME OR FORTRAN COMPILER IN VAX SYSTEM 
C 


IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,KSIC,MUO 

DIMENSION X(IOO) ,Y(100) .ERROR (100) .FEEREC(IOO) .UREC(IOO) , 
+ THETAR (100) 


DEFINE PARAMETERS USED BY PROGRAMS 

(1) IN ANG LP ARE UNIT NUMBERS ASSIGNED TO THE INPUT AND OUTPUT 
DEVICES 

IN=5 
LP=6 

C (2) NDBUG IS USE TO CONTROL THE AUXILIARY OUTPUT FOR DEBUGGING 
NDBUG=1 

C (3) NSOLVE IS THE UPPER LIMITATION OF REPEATATION FOR SOLVING 
C NONLINEAR EQUTIONS; 

C EPSI IS THE CLEARANCE OF FUNCTION VALUES WHEN THE FUNCTIONS 

C IS CONSIDERED AS SOLVED (ALL FUNTIONS HAVE FORMS OF F(X)=0); 

C DELTA IS THE CLERANCE OF VARIABLE INCREMENT WHEN FUNCTION IS 

C SOLVED 

C NC, EPSI AND DELTA MAY BE CHANGED WHEN SOLUTIONS ARE DIVERGENT 

C OR LESS ACCURATE 

NSOLVE=100 
DELTA=1.D-15 
EPSI=1.D-15 

C (4) OTHER PARAMETERS (DON'T CHANGE) 

DR=DATAN ( 1 . DO) /45 . DO 
C 

C DEFINE INPUT PARAMTERS OF PROBLEM (USE INCH AS UNIT OF LENGTH) 

C 

C (1) PINION AND GEAR: PN=DIAMETRAL PITCH; Nl=PINION TOOTH NUMBER; 

C KS IN=PRESSURE ANGLE IN NORMAL SECTION (DEGREE) ; 

C BETAP=HELIX ANGLE (DEGREE) ; 

C HD=HEIGHT OF DEDENDUM OF PINION; 
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noon non no 


HA=HEIRHT OF ADDENDUM OF PINION 

ZCOE=COEFFICIENT OF PINION TOOTH LENGTH (THE LENGTH" ZCOE/PN) 
PN=10.D0 
Nl=20 

KS IN=20 . D0*DR 
BETAP=20.D0*DR 
HD=1 .DO/PN 
HA=1 .DO/PN 
ZCOE=10. 

(2) TOOL: ALPHA=HALF OF CONE VERTEX ANGLE (DEGREE) ; 

RC=RADIUS OF BOTTOM CIRCLE OF CONE; 

MU=TILT ANGLE TO INSTALL PINION-CUTTING TOOL 

ALP=80.D0*DR 

MU0=DATAN (DSIN (KS IN) *DTAN (BETAP) ) 

MU=1.*MU0 
RC= 1 . DO 

(3) OUTPUT: 

NZ=NO. OF CROSS SECTIONS WHERE PINION PROFILE IS SIMULATED; 
NU=NO. OF MAXIMUM POINTS USED FOR SIMULATE PINION PROFILE 
UIN= INCREMENT OF PINION TOOTH SURFACE COORDINATE U 
NZ=21 
NU=61 

UIN=2 . 5D0*CH/ (NU- 1) 

C 

C DESCRIPTION OF OUTPUT VARIABLES 

C Z1=DISTANCE BETWEEN CROSS SECTION CONSIDERED AND MIDDLE CROSS 
C SECTION 

C NO=OUTPUT NO. 

C U=TOOL SURFACE COORDINATE 

C THETA=TOOL SURFACE COORDINATE 

C FEE'PINION TOOTH SURFACE GENERATION PARAMETER 

C X1=X COORDINATE OF PINION PROFILE 

C Y1=Y COORDINATE OF PINION PROFILE 

C Rl'RADIUS OF PINION PROFILE 

C VSH= AVERAGE DEVIATION SHIFT OF CROSS SECTION PROFILE FROM PROFILE 

C OF GENERAL PINION (INVOLUTE CURVE) 

C VPE=MAXIMUM DEVIATION OF CROSS SECTION PROFILE FROM INVOLUTE CURVE 
C VSD= STANDARD DEVIATION OF CROSS SECTION PROFILE FROM INVOLUTE 

C CURVE 

C 

C THE PROGRAM IS WRITTEN BY JIAO ZHANG 
SK=DSIN (KSIN) 

CK=DCOS (KSIN) 

SB=DS IN (BETAP) 

CB=DCOS (BETAP) 

SM=DSIN (MU) 

CM"DCOS (MU) 

SA=DS IN (ALPHA) 

CA=DCOS (ALPHA) 

KSIC=DATAN (SK/CK/CB) 

CKC=DCOS (KSIC) 

SKC=DSIN (KSIC) 

PT=PN*CB 

RP=Nl/2./PT 
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RB=RP*DCOS (KSIC) 

RA=RP+HA 

CL=RC/SA 

CH=HD/CK/CM 

Al=CA*CK*CB+SA*SK*CM*CB+SA*SM*SB 

A2=SK*SM*CB-CM*SB 

A3=CA*SK*CM*CB-SA*CK*CB+CA*SM*SB 

A4=SK*CM*CB+SM*SB 

B1=CA*SK-SA*CK*CM 

B2=CK*SM 

B3=SA*SK+CA*CK*CM 

B4=CK*CM 

C 1 =CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 

C2=SK*SM*SB+CM*CB 

C3=CA*SK*CM*SB-SA*CK*SB~CA*SM*CB 

C4=-SK*CM*SB+SM*CB 

D1=CM*CB+SK*SM*SB 

D2=SA*SM*CB-SB* (CA*CK+SA*SK*CM) 

D3=CA*CK*SB 

D4=CA* (CA*SK-SA*CK*CM) 

D5=SA* (SA*SK+CA*CK*CM) 

D6=CA*CK*SM 
DO 5 1=1, NZ 

Zl=-ZCOE/PN/2 . +ZCOE/PN*FLOAT (1-1) /FLOAT (NZ-1) 

FEEO=Zl*SB/CB/RP 
CFO=DCOS (FEEO) 

SFO=DSIN (FEEO) 

R1=RP 

ERR=0. 

NP=0 

DO 15 J= 1 , NU 

IF (Rl.GT.RA.AND. J.GT.l) GOTO 55 
U=CL-UIN*FLOAT ( J-l) 

TEMP1= (Zl- (CL-CH) *C4-U*CA*C3) /U/SA/DSQRT (C1*C1+C2*C2) 
TEMP2=DARS IN (C 1 /DSQRT (C 1 *C 1 +C2*C2) ) 

TEMP 3=DARS IN (TEMP 1 ) 

THE TA= TEMP 3- TEMP 2 
CT=DCOS (THETA) 

ST=DS IN (THETA) 

XC=U*SA* (CT*A1+ST*A2) +U*CA*A3~ (CL-CH) *A4 
YC=U*SA* (CT*B1-ST*B2) -U*CA*B3+ (CL-CH) *B4 

FEE= (U* (CT*D1+ST*D2) - (CL-CH) * ( (CT*CA*CA+SA*SA) *D1~ST*D3) ) /RP/ 
+ (CT*D4+D5-ST*D6) 

CF=DCOS (FEE) 

SF=DSIN (FEE) 

X1=XC*CF+YC*SF-RP*FEE*CF+RP*SF 
Y1=-XC*SF+YC*CF+RP*FEE*SF+RP*CF 
R1=DSQRT (X1*X1+Y1*Y1) 

IF (Rl.LT.RB) GOTO 15 

NP=NP+1 

X(NP)=X1 

Y(NP)=Y1 

FEEREC (NP) =FEE 

UREC (NP) =U 
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n n 


THETAR (NP) “THETA 
XE=X 1 *CFO+ Y 1 *SFO 
YE=-X1*SF0+Y1*CF0 
YS=YE 

FEET=FEE+FEEO 
DO 25 K=l,NSOLVE 

W=RP*FEET*CKC*DS IN (FEET-KS IC) +RP*DCOS (FEET) - YS 
DWDF=-RP*SKC*DCOS (FEET-KSIC) +RP*FEET*CK*DCOS (FEET-KSIC) 
DF=-W/DWDF 

IF (DABS (W) .LT. EPS I. AND. DABS (DF) .LT. DELTA) GOTO 65 
FEET=FEET+DF 
25 CONTINUE 
65 CONTINUE 

XS=-RP*FEET*CKC*DCOS (FEET-KSIC) +RP*DSIN (FEET) 

ERROR (NP) =XE-XS 

WRITE (LP.110) K, W,DF, ERROR (NP) , FEET, FEE+FEEO 
110 FORMAT (IX, ' K= ' ,I5,5X, ' W= ' , D15 . 7, 5X, ' DF= ' ,D15 . 7,D15 . 7 , 2D25 . 17) 
ERR=ERR+ERROR (NP) 

15 CONTINUE 
55 VSH=-ERR/FLOAT(NP) 

VPE=0. 

VSD=0. 

DO 35 J=1 , NP 
ETEMP=ERROR ( J) +VSH 

IF (DABS(ETEMP) .GT.VPE) VPE=DABS (ETEMP) 

VSD=VSD+ETEMP**2 
35 CONTINUE 

VSD=DSQRT (VSD/FLOAT (NP) ) 

WRITE (LP, 10) Z1 , VSH, VPE.VSD 

10 FORMAT (1H1,///'Z1=' ,F15.7, 5X, ' VSH= ' , F15 . 7 , 5X, ' VPE- ' ,F15. 7,5X, 

+ ' VSD ' ,D15.7,5X/2X, 'NO. ' ,7X, ' U ' , 14X, ' THETA' , 10X, ' FEE ' , 12X, 

+ 'XI' ,13X, ' Y1 ' , 13X, ' R1 ' , 13X, 'ERROR') 

DO 45 J=1,NP 

Rl=DSQRT (X ( J) **2+Y ( J) **2) 

ETEMP=ERROR ( J) +VSH 

WRITE (LP , 20) J,UREC (J) , THETAR ( J) , FEEREC (J) ,X (J) , Y (j) ,Rl , ETEMP 
20 FORMAT (2X, 12, 2X, 7F15 . 7) 

45 CONTINUE 
5 CONTINUE 
STOP 
END 

SENTRY 

C FIND AUXILIARY VALUES FOR CALCULATION 
RP=FLOAT (Nl) /2.D0/PN 
CL=RC/DSIN (alp) 

D=CL*DCOS (ALP) 

Al=CL-HDC/DCOS (RKS) 

RPU=RP+HAC 
DO 5 1=1, NL 
Z=FL0AT(I-1)*ZI 
YY=D-Z/DTAN(ALP) 

WRITE (LP, 10) Z 

10 FORMAT (1H1/1X, ' Zl= ' .F15.7/1X, 'NO' , 10X, ' Yl ' , 13X, 'XP' ,13X, ' YP ’ , 

+ 13X, 'R1 ’ ,13X, 'FEE') 
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C CALCULATE THE PROPILE OF THE SURFACE CUT BY PLANE Z=CONST 
KKK=0 

DO 15 J=1,N 

Yl=YY*FLOAT ( J-l) /FLOAT (N~l) 

F=DSQRT (1 . - (Z/ (D-Yl) /DTAN (ALP) ) **2) 

FEE= ( (D-Yl) *F/DCOS (ALP) -Al* (F*DCOS (ALP) *DCOS (ALP) +DSIN (ALP) *DSIN ( 
+ ALP) )) /RP/ (DSIN (ALP) *DCOS (ALP-RKS) -F*DCOS (ALP) *DSIN (ALP-RKS) ) 
XP (I , J) = (D-Yl) * (DTAN (ALP) *F*DCOS (ALP-RKS+FEE) -DSIN (ALP-RKS+FEE) ) 

+ +A1 *DS IN (FEE-RKS) -RP*FEE*DCOS (FEE) +RP*DSIN (FEE) 

YP (I , J) =- (D-Yl) * (DTAN (ALP) *F*DSIN (ALP-RKS+FEE) +DCOS (ALP-RKS+FEE) ) 
+ +Al*DCOS (FEE-RKS) +RP*FEE*DSIN (FEE) +RP*DCOS (FEE) 

Rl=DSQRT (XP (I , J) **2+YP (I , J) **2) 

IF (Rl.GT.RPU) THEN 
JJ=J-1 

XP (I , J) =XP (I , J J) + (XP (I , J) -XP (I , J J) ) / (R1-R1TEMP) * (RPU-R1TEMP) 

YP (I , J) =YP (I , J J) + (YP (I , J)-YP (I , J J) ) / (R1-R1TEMP) * (RPU-R1TEMP) 
FEE=FEETEM+ (FEE-FEETEM) / (R1-R1TEMP) * (RPU-R1TEHP) 

Y1=Y1TEM+ (Y1-Y1TEM) / (R1-R1TEMP) * (RPU-RlTEMP) 

Rl=DSQRT (XP (I , J) **2+YP (I , J) **2 ) 

KKK=1 
END IF 

WRITE (LP , 20) J,Y1,XP(I,J) , YP (I , J) ,R1, FEE 
20 FORMAT (1X,I4,^F15.7,F15.7) 

IF (KKK.GT.O) GO TO 30 
FEETEM=FEE 
Y1TEM=Y1 
R1TEMP=R1 
15 CONTINUE 
30 NS (I) = J-l 

IF (I.NE.1) GO TO 55 
NS1=NS(1) 

GO TO 5 

C PREPARATION OF INTERPLORATION 
55 NS2=NS(I) 

DO 105 L=1 ,NS2 
J«NS2+1-L 

IF (YP(1,NS1) .GT.YP(I, J)) GO TO 110 
105 CONTINUE 
110 NS2-J 
NREC=2 
XERS=O.DO 
DO 115 L-1.NS2 
DO 125 J=NREC,NS1 
IF (YP(1, J) .GT.YP(I.L)) GO TO 120 
125 CONTINUE 
120 NREC'J 
J1=J-1 

XERROR(L) =XP(I,L)-XP(1,J1)-(YP(I,L)-YP(1,J1))*(XP(1,J)-XP(1,J1)) 

+ / (YP (1 , J) -YP (1 , Jl) ) 

XERS=XERS+XERROR (L) 

115 CONTINUE 

XERS=-XERS/FLOAT (NS2) 

XPE=0.D0 

SDX=0.D0 


99 



IF (NDBUG.GT.2) WRITE (LP,80) 

80 FORMAT (1H1 , 13X, 'NO. ' , 4X, 'DIVIATION VALUE') 

DO 45 L=1 ,NS2 

XERROR (L) =XERROR (L) +XERS 

IF (NDBUG.GT.2) WRITE (LP,50) L, XERROR (L) 

50 FORMAT (13X, I3.2F15.7) 

IF (DABS (XERROR (L) ) . GT . XPE) XPE=DABS (XERROR (L) ) 
SDX=SDX+XERROR (L) **2 
45 CONTINUE 

SDX=DSQRT (SDX/FLOAT (NS2) ) 

WRITE (LP.40) XERS ,XPE, SDX 

40 FORMAT (///IX, 'XSH* ' , E15 . 7,5X, 'XPE® ' , E15. 7,5X, ' SDX*= ' ,E15. 7) 
IF (NDBUG.GT.2) WRITE (LP,60) NS1.NS2 
60 FORMAT (//5X, ' NS1= ' , 16 , 6X, ' NS2= ' , 16) 

5 CONTINUE 
STOP 
END 
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C. 

C. 

C. 

C. 

C. 

C. 

C. 

C. 

C. 

C. 


******************************************************* 


* 


* 


* PROGRAM I I 

* UDDERCUTTING CONDITION FOR HELICAL PINION 

* GENERATED BY CONE CUTTER 

* AUTHOR: FAYDOR LITVIN 

* JIAO ZHANG 

* 

* 


* 

* 

* 

* 

* 

* 

* 


******************************************************* 


C 

C PURPOSE 
C 

C THIS PROGRAM IS USED TO FIND THE UDDERCUTTING CONDITIONS FOR A 
C HELICAL PINION GENERATED BY CONE CUTTER 
C 

C NOTE 
C 

C THIS PROGRAM IS WRITTEN IN FORTRAN 77. IT CAN BE COMPILED BY V 
C COMPILE IN IBM MAINFRAME OR FORTRAN COMPILER IN VAX SYSTEM. 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN.MU 


DEFINE PARAMETERS USED BY PROGRAMS 


(1) IN ANG LP ARE UNIT NUMBERS ASSIGNED TO THE INPUT AND OUTPUT 

DEVICES 

IN=5 

LP=6 

(2) NDBUG IS USE TO CONTROL THE AUXILIARY OUTPUT FOR DEBUGGING 
NDBUG=2 

(3) OTHER PARAMETERS (DON'T CHANGE) 

DR=DATAN ( 1 . DO) / 45 . DO 

DEFINE INPUT PARAMTERS OF PROBLEM (USE INCH AS UNIT OF LENGTH) 

(1) PINION AND GEAR: PN=DIAMETRAL PITCH; Nl=PINION TOOTH NUMBER; 

KSIN=PRESSURE ANGLE; BETAP=HELIX ANGEL OF PINION; 

HD=HEIGHT OF DEDENDUM OF PINION 
PN-10.D0 
Nl=20 

KSIN=20.D0*DR 
BETAP=30.D0*DR 
HD=1 .DO/PN 

(2) TOOL: ALPHA=HALF OF CONE VERTEX ANGLE (DEGREE) *, 

RC=RADIUS OF BOTTOM CIRCLE OF CONE; 

MU=TILT ANGLE OF MOUNTING TOOL 
ALPHA=20.D0*DR 
RC=0. 35D0 

MUO-DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU=0 . *MUO 

(3) PROBLEM: NPROB=ID NO. OF PROBLEM (-1=GIVEN Nl AND HD, FIND IF 

UDDERCUTTING OCCUR; 0=GIVEN Nl, FIND MAXIMUM HD WITHOUT UDDER- 
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CUTTING; 1=GIVEN HD, FIND MINIMUM Nl WITHOUT UDDERCUTTING); 
N=NO. OF THETA VALUES USED CALCULATION; 

DELTHE= INCREMENT OF THETA (DEGREE) 

NPROB= 1 
N=21 

DELTHE=l.DO*DR 
C 

C DESCRIPTION OF OUTPUT 
C 

C OUTPUT IS A STATEMENT BASED ON THE PROBLEM WITHOUT ANY LITERAL 
C PARAMETER 
C 

C FIND AUXILIARY VALUES FOR CALCULATION 
SK=DS IN (KS IN) 

CK=DCOS (KSIN) 

SB=DSIN (BETAP) 

CB=DCOS (BETAP) 

SM=DSIN (MU) 

CM=DCOS (MU) 

SA=DS IN (ALPHA) 

CA=DCOS (ALPHA) 

PT=PN*CB 

RP=Nl/2./PT 

CL=RC/SA 

CH=HD/CK/CM 

UU=CL~CH 

AA=CM*CB+SK*SM*SB 

BB=CA*CK*SB+SA*SK*CM*SB~SA*SM*CB 

CC=CA*CK*SB 

DD1=CA*SK~SA*CK*CM 

DD=DD1*CA 

EE 1 = SA*SK+CA*CK*CM 

EE=EE1*SA 

FF1=CK*SM 

FF=FF1*CA 

11= (CA*SK*CM*SB-SA*CK*SB-CA*SM*CB) *CA/SA 
WRITE (LP,90) 

90 FORMAT (1H1) 

IF (NPROB) 5,15,25 
C CHECK IF UNDERCUTTING OCCURS 
5 WRITE (LP, 10) 

10 FORMAT (1H1/3X, 'NO' ,7X, ’THETA' ,12X, 'A' ,13X, 'B' ,14X, 'C' ,10X, 

+ 'b**2-4.*A*C' ,3X, 'U1/(RC/S IN (ALPHA)) ') 

UMIN=5.D0 
DO 45 1=1, N 
NN= (N+l)/2 

THE=DBLE (FLOAT (I-NN) ) *DELTHE 
ST=DSIN (THE) 

CT=DCOS (THE) 

GG=DD*CT+EE-FF*ST 

WW= (BB*BB+AA*AA) *EE+ (BB*DD-AA*FF) *1 I+ST*ST*ST* (AA*AA*FF-BB*BB*FF 
+ + 2 . *AA*BB*DD) -CT*CT*CT* (BB*BB*DD-AA*AA*DD+2 . *AA*BB*FF) 

+ -ST* (2 . *AA*AA*FF+AA*BB*DD-AA*EE*II) +CT* (2 . *BB*BB*DD+AA*BB*FF 

+ +BB*EE*II) 
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XX= ( (AA*FF*SA*SA-CC*EE) *CT+ (AA*DD*SA*SA-AA*EE*CA*CA) *ST+ (AA*FF*CA 
+ *CA-CC*DD))*(BB*CT+AA*ST+Il) 

YY=DD1 *SA*CT-FF 1 *SA*ST~EE 1 *CA 

ZZ=CM*CK 

A=YY*WW 

B=UU* (WW*ZZ+XX*YY) +RP*GG*GG*GG 

C=UU*UU*XX*ZZ 

d=b*b-4.*a*c 

IF (D.GE.O.) THEN 

Ul= (-B-DSQRT (D) ) /2. /A 

U2= (-B+DSQRT (D) )/2./k 

U1=U1/CL 

IF (UMIN.GT.U1) UMIN=U1 

U2=U2/CL 

THE=THE/DR 

IF (NDBUG.LT. 1) WRITE (LP.lOO) I,THE,A,B,C,D,Ul 
100 FORMAT (IX, I4,8F15. 7) 

ELSE 

IF (NDBUG.LT. 1) WRITE (LP,100) I , THE, A, B , C, DD 
END IF 
45 CONTINUE 

IF (UMIN.LT.1.D0) WRITE (LP,110) 

110 FORMAT (///IX, 'UDDERCUTTING WILL OCCUR FOR YOUR DESIGN') 

IF (UMIN.GE. 1 .DO) WRITE (LP.120) 

120 FORMAT (///IX, ' UDDERCUTTING WILL NOT OCCUR FOR YOUR DESIGN') 

GO TO 35 

C DETERMINE THE MAXIMUM ADDENDUM HEIGHT OF RACK CUTTER 
15 WRITE (LP , 20) 

20 FORMAT (1H1/3X, ' NO ' , 7X, ' THETA' , 12X, ' A ' , 13X, ' B ' , 14X, ' C ' , 10X, 

+ 'B**2-4.*A*C' ,3X, 'ALLOWED RATIO OF HD/(1/PN)') 

UMIN=5 . DO 
DO 55 1=1, N 

THE=DBLE (FLOAT (I-NN) ) *DELTHE 
ST=DSIN (THE) 

CT=DCOS (THE) 

GG=DD*CT+EE-FF*ST 

WW= (BB*BB+AA*AA) *EE+ (BB*DD~AA*FF) *II+ST*ST*ST* (AA*AA*FF-BB*BB*FF 
+ +2 . *AA*BB*DD) -CT*CT*CT* (BB*BB*DD-AA*AA*DD+2 . *AA*BB*FF) 

+ -ST* (2 . *AA*AA*FF+AA*BB*DD~AA*EE*I I) +CT* (2 . *BB*BB*DD+AA*BB*FF 

+ +BB*EE*II) 

XX= ( (AA*FF*SA*SA-CC*EE) *CT+ (AA*DD*SA*SA-AA*EE*CA*CA) *ST+ (AA*FF*CA 
+ *CA-CC*DD) ) * (BB*CT+AA*ST+II) 

YY=DD1 *SA*CT-FF 1 *SA*ST~EE 1 *CA 

ZZ=CM*CK 

A=XX*ZZ 

B=CL* (WW*ZZ+XX*YY) 

C=CL*CL*YY*WW+CL*RP*GG*GG*GG 

D=B*B-4.*A*C 

TEST=DABS (A/B) 

EPS=1.D-16 
THE= THE/DR 
IF (D.GE.O.) THEN 
IF (TEST.GT.EPS) THEN 
Ul= (-B+DSQRT (D) ) /2. /A 
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U2= (-B-DSQRT (D) ) /2 . /A 

ELSE 

Ul=-C/B 

U2=0. 

END IF 

Ul= (CL-U1) *CK*PN 
U2= (CL-U2) *CK*PN 

IF (NDBUG.LT. 1) WRITE (LP.100) I,THE,A,B,C,D,U1,U2 
ELSE 

IF (NDBUG.LT. 1) WRITE (LP.lOO) I,THE,A,B,C,D 
END IF 

IF (UMIN.GT.U1) UMIN=U1 
55 CONTINUE 

WRITE (LP , 200) UMIN 

200 FORMAT (///IX, 'TO AVOID UDDERCUTTING, IT IS NECESSARY TO KEEP DEDE 
+NDUM OF PINION <= ' , F10. 7 , ' /PN ' ) 

GO TO 35 

C DETERMINE THE MINIMUM NO. OF TEETH FOR UNUNDERCUTTING 
25 WRITE (LP , 30) 

30 FORMAT (1H1/3X, 'NO' ,7X, 'THETA' ,7X, 'NO. OF TEETH') 

RNMAX=0. 

DO 65 1=1, N 

THE=DBLE (FLOAT (I-NN) ) *DELTHE 
ST=DSIN (THE) 

CT=DCOS (THE) 

GG=DD*CT+EE-FF*ST 

WW= (BB*BB+AA*AA) *EE+ (BB*DD“AA*FF) *II+ST*ST*ST* (AA*AA*FF-BB*BB*FF 
+ +2 . *AA*BB*DD) -CT*CT*CT* ( B B *B B *DD~AA* AA*DD+ 2 . *AA*BB*FF) 

+ -ST* (2 . *AA*AA*FF+AA*BB*DD-AA*EE*II) +CT* (2 . *BB*BB*DD+AA*BB*FF 

+ +BB*EE*II) 

XX= ( (AA*FF*SA*SA-CC*EE) *CT+ (AA*DD*SA*SA-AA*EE*CA*CA) *ST+ (AA*FF*CA 
+ *CA-CC*DD))*(BB*CT+AA*ST+II) 

Y Y=DD1 *SA*CT-FF 1 *SA*ST~EE 1 *CA 
ZZ=CM*CK 

RR=- (CL*CL*YY*WW+CL*UU* (WW*ZZ+XX*YY) +UU*UU*XX*ZZ) /CL/GG**3 

RN=2.*RR*PT 

THE=THE/DR 

IF (RNMAX.LT. RN) RNMAX=RN 
IF (NDBUG.LT. 1) WRITE (LP,'100) I , THE , RN 
65 CONTINUE 

WRITE (LP , 300) RNMAX 

300 FORMAT (/// IX, 'WITHOUT UDDERCUTTING, MINIMUM TOOTH NO. OF PINION 
+IS: ' ,F11.7) 

35 WRITE (LP , 400) 

400 FORMAT (1H1) 

STOP 

END 
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******************************************************* 


C. . 
C. . 

c. . 
c. . 
c. . 
c. . 
c. . 


c. . 
c. . 
c.. 
c. . 


* * 

* PROGRAM III * 

* CONTACT ELLIPSIS FOR HELICAL PINION GENERATED BY * 

* CONE CUTTER IN MESHING WITH REGULAR HELICAL GEAR * 

* * 

* AUTHORS: FAYDOR LITVIN * 

* JIAO ZHANG * 

* * 

* * 


********************************** ********************* 


C 

C PURPOSE 
C 

C THIS PROGRAM IS USED TO FIND THE SHAPE AND ORIENTAION OF THE CONTACT 
C ELLIPSE WHEN A HELICAL PINION CROWNED BY CONE CUTTER IS IN 

C MESHING WITH A REGULAR HELICAL GEAR 

C 

C NOTE 
C 

C THIS PROGRAM IS WRITTEN IN FORTRAN 77. IT CAN BE COMPILED BY V 
C COMPILER IN IBM MAINFRAME OR FORTRAN COMPILER IN VAX SYSTEM. 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,MUO,MUT,KSIC,KF,KH,KSP,KQP,KSG,KQG 
DIMENSION CMM(3,4) ,EFD(3) ,EHD(3) ,RND(3) ,RD(4) ,EFF(3) ,EHF(3) , 

+ RNF(3) ,RC(3) ,RF(3) ,Rl (3) 


C 

C DEFINE PARAMETERS USED BY PROGRAMS 
C 

C (1) IN ANG LP ARE UNIT NUMBERS ASSIGNED TO THE INlUT AND OUTPUT 
C DEVICES 


IN=5 


LP=6 

C (2) NDBUG IS USE TO CONTROL THE AUXILIARY OUTPUT FOR DEBUGGING 
NDBUG=2 

C (3) OTHER PARAMETERS (DON ' T CHANGE) 

DR=DATAN ( 1 . DO) /45 . DO 
C 

C DEFINE INPUT PARAMTERS OF PROBLEM (USE INCH AS UNIT OF LENGTH) 

C 

C (1) PINION AND GEAR: PN'DIAMETRAL PITCH; Nl=PINION TOOTH NUMBER; 
C MPG=TOOTH NUMBER RATIO (GEAR TOOTH NO./N1); 

C KS IN=PRESSURE ANGLE IN NORMAL SECTION; 

C BETAP=HELIX ANGEL; 

C HD=HEIGHT OF DEDENDUM OF PINION 

PN=10.D0 
Nl=20 
MPG=2 

KSIN=20.D0*DR 
BETAP*= 1 0 . D0*DR 
HD=1./PN 

(2) TOOL: ALPHA“HALF OF CONE VERTEX ANGLE (DEGREE) ; 

RC=RADIUS OF BOTTOM CIRCLE OF CONE; 
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C MU=TILT ANGEL TO INSTALL PINION CUTTING TOOL 

ALPHA=80.D0*DR 

MUO=DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU=0.*MU0 
RC=1 .DO 

C (3) DEFORMATION: DEL=CONTACT DEFORMATION AT CONTACT POINT 
DEL=4.D-4 

(4) OUTPUT: NU=NUMBER OF CONTACT POINTS IN MATING SURFACES FOR US 
TO CALCULATE CONTACT ELLIPSES) 

NU=101 
C 

C DESCRIPTION OF OUTPUT PARAMETER 
C 

C Rl=PINION RADIUS OF CONTACT POINT 

C ALPHA=THE ROTATION ANGLE BETWEEN PRINCIPAL DIRECTION OF PINION 
C TOOTH SURFACE AND AXES OF CONTACT ELLIPSE 

C A=LENGTH OF HALF SHORT AXIS OF CONTACT ELLIPSE 

C B=LENTH OF HALF LONG AXIS OF CONTACT ELLIPSE (ALONG DIRECTION OF 
C GEAR TOOTH LEHGTH) 

C RNF=UNIT NORMAL OF PINION TOOTH SURFACE AT CONTACT POINT 

C EFF=PRINCIPAL DIRECTION OF PINION TOOTH SURFACE AT CONTACT POINT 

C EHF=PRINCIPAL DIRECTION OF PINION TOOTH SURFACE AT CONTACT POINT 

C 

C FIND AUXILIARY VALUES FOR CALCULATION 
SK=DSIN (KSIN) 

CK=DCOS (KSIN) 

SB=DS IN (BETAP) 

CB=DCOS (BETAP) 

SM=DSIN (MU) 

CM=DCOS (MU) 

SA=DS IN (ALPHA) 

CA=DCOS (ALPHA) 

KSIC=DATAN (SK/CK/CB) 

CKC=DCOS (KSIC) 

SKC=DSIN (KSIC) 

PT=PN*CB 
RP=Nl/2./PT 
RB=RP*DCOS (KSIC) 

RA=RP+1./PN 

N2=N1*MPG 

RG=RP*MPG 

CL=RC/SA 

CH=HD/CK/CM 

A=CL-CH 

CMM(1 , 1) =CA*CK*CB+SA*SK*CM*CB+SA*SM*SB 
CMM (1 , 2) =SA*CK*CB-CA*SK*CM*CB-CA*SM*SB 
CMM (1 , 3) =SK*SM*CB-CM*SB 
CMM (1,4)=- (SK*CM*CB+SM*SB) 

CMM (2 , 1) =CA*SK-SA*CK*CM 
CMM (2 , 2) =SA*SK+CA*CK*CM 
CMM (2 , 3) =-CK*SM 
CMM (2 , 4) =CK*CM 

CMM (3 , 1) =CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 
CMM (3 , 2) =SA*CK*SB-CA*SK*CM*SB+CA*SM*CB 
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CMM (3 , 3) =SK*SM*SB+CM*CB 
CMM (3 , 4) =-SK*CM*SB+SM*CB 
D1=CM*CB+SK*SM*SB 
D2=SA*SM*CB~SB* (CA*CK+SA*SK*CM) 

D3=CA*CK*SB 

D4=CA* (CA*SK-SA*CK*CM) 

D5=SA* (SA*SK+CA*CK*CM) 

D6=CA*CK*SM 

UL=CL-2.D0*CH 

UU=CL 

MUT=MUO/DR 

WRITE (LP ,110) RP,RA,RB,MUT 

110 FORMAT (///IX, ’ RP= 1 , F15 . 7, 5X, 1 RA= ’ , F15 . 7 ,5X, 1 RB= 1 , F15.7, 5X, 

+ 1 MUO= 1 , F15 . 7) 

TH=0. DO 
ST=DSIN(TH*DR) 

CT=DCOS (TH*DR) 

EFD (1) =ST 

EFD(2)=0.D0 

EFD(3)=-CT 

EHD(1)=-CT*SA 

EHD(2)=CA 

EHD(3)=-ST*SA 

RND (1) =CT*CA 

RND(2)=SA 

RND (3) =ST*CA 

CALL MATMUL (CMM, EFD, EFF, 3 , 3) 

CALL MATMUL (CMM.EHD, EHF, 3 , 3) 

CALL MATMUL (CMM, RND,RNF , 3 , 3) 

WRITE (LP, 10) TH, (EFF(M) ,M=1,3) , (EHF(M) ,M=1,3) , (RNF(M) ,M=1,3) 

10 FORMAT (1H1.///3X, ' THTEA: ' ,F15.7, ' DEGREE '//IX, ' EFF: ’ ,3F15.7/1X, 

+ 'EHF: ' .3F15.7/1X, 'NF :',3F15.7) 

DO 15 J=1 , NU 

U=UU- (UU-UL) /FLOAT (NU~ 1) *FLOAT ( J- 1) 

KF=-CA/SA/U 

KH=O.DO 

RD(1)=U*CT*SA 

RD(2)=-U*CA 

RD(3) =U*ST*SA 

RD(4)=A 

CALL MATMUL (CMM , RD , RC , 3 , 4) 

FEE= (U* (CT*D1+ST*D2) -A* ( (CT*CA*CA+SA*SA> *D1~ST*D3) ) / (CT*D4+D5-ST* 
+ D6) /RP 

RF (1) =RC (1) -RP*FEE 
RF (2) =RC (2) +RP 
RF (3) =RC (3) 

RPP=DSQRT (RF (1) **2+RF (2) **2) 

CALL PR0T(RF,R1,FEE) 

TEMP=RC (2) *EFF (1) -RF (1) *EFF (2) 

B13=EHF(3)-KF*TEMP 

B23=~EFF(3) 

B33=- (RNF (1) *RF (1) +RNF (2) *RF (2) +KF*TEMP*TEMP) 

SIGMA=0. 5D0*DATAN (2 . *B 13*B23/ (B23*B23~B 13*B 13~ (KF-KH) *B33) ) 

C0E1= (B23*B23-B13*B13- (KF-KH) *B33) /B33/DCOS (2. *SIGMA) 
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C0E2=KF+KH+ (B13*B13+B23*B23) /B33 
KSP= (C0E2-C0E1) /2 . DO 
KQP= (C0E2+C0E1) /2.D0 

IF (NDBUG.LT. 1) WRITE (LP,20) U,FEE,RPP, (R1 (M) ,M=1,3) , 

+ B13,B23,B33, SIGMA, KSP,KQP 

20 FORMAT (//IX, 'U = ' , F15 . 7,5X, ' FEE= ' , F15 . 7 , 5X, 'R1 =',F15.7/ 

+ IX, 'XP =' ,F15.7,5X, ' YP =' ,F15.7,5X, 'ZP = ' , F15 . 7/ 

+ IX, 'B13=' ,F15.7,5X, ' B23= ' , F15 . 7 ,5X, ' B33= ' , F15 . 7/ 

+ IX, 1 SIG= ' ,F15.7,5X, 1 KSP= ' ,F15.7,5X, 'KQP=' ,F15.7) 

W=- (A-U) 

KSG=0. 

KQG=1 . / (RG*SK* (CKC/CK/CB) **2+W*CM*CK/SK) 

G1=KSP-KQP 

G2=KSG-KQG 

S1=KSP+KQP 

S2=KSG+KQG 

SIGMGP=-SIGMA 

ALPHA1=0. 5D0*DATAN (G2*DSIN (2. *SIGMGP) / (G1-G2*DC0S (2. *SIGMGP) ) ) 
ALPHA2=ALPHA1+SIGMGP 

AA= (S1-S2-DSQRT (Gl*Gl-2. *G1*G2*DC0S (2 . *SIGMGP) +G2*G2) ) /A . 

BB= (S1-S2+DSQRT (Gl*Gl-2. *G1*G2*DC0S (2. *SIGMGP) +G2*G2) ) /A. 

AAA=1 . /DSQRT (DABS (AA) ) 

BBB=1 . /DSQRT (DABS (BB) ) 
ratio=bbb/aaa 

ALPHAl=ALPHAl/DR 

ALPHA2=ALPHA2/DR 

WRITE (LP , 130) G1 , G2 , S 1 , S2 , ALPHAl 
130 FORMAT (1X.5F15.7) 

WRITE (LP , 30) RPP,ALPHA2, AAA, BBB, RATIO 
30 FORMAT (IX, 'Rl = ' , F15 . 7 ,5X, ' ALP= 1 , F15 . 7, 5X, ' A =',F15.7,5X, 

+ 'B =' ,F15.7,5X, 'B/A=' ,F15.7) 

IF (RPP.GT.RA) GO TO 5 
15 CONTINUE 
5 STOP 
END 


SUBROUTINE MATMUL(CMM,A,B,N,M) 

THIS SUBROUTINE IS USED TO MULTIPLY THE MATRIX CMM(N*M) BY THE MATRIX 
A(M*1) . THE RESULT IS STORED IN THE MATRIX B(N*1) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION CMM(3,M) ,A(M) ,B(N) 

DO 5 1=1, N 
5 B (I) =0. 

DO 15 1 = 1, N 
DO 15 J=1,M 

15 B(I)=B(I)+CMM(I, J)*A(J) 

RETURN 

END 


SUBROUT INE PROT (A , B , FEE) 

THIS SUBROUTINE IS USED TO ROTATE COORDINATE SYSTEM PLANARLY IN XOY 
THROUGH ANGLE FEE 
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DOUBLE PRECISION A (3) ,B(3) ,FEE 
B (1) =A(1) *DCOS (FEE) +A(2) *DSIN (FEE) 
B (2) =A (2) *DCOS (FEE) -A ( 1 ) *DS IN (FEE) 
B(3)=A(3) 

RETURN 

END 
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PROGRAM IV 

TRANSMISSION ERRORS OF HELICAL PINION GENERATED 
BY CUTTER WITH CONE OR REVOLUTE SURFACE 
IN MESHING WITH REGULAR HELICAL GEAR 


AUTHORS: 


FAYDOR LITVIN 
JIAO ZHANG 


* 

* 

* 

* 

A 

A 

A 

A 

A 

A 

A 


AAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAA AAAAA AAAAAAAAA A 


c 

C PURPOSE 
C 

C THIS PROGRAM IS USED TO CALCULATE THE TRANSMISSION ERRORS OF A 
C PINION GENERATED BY THE CUTTER WITH CONE OR REVOLUTE SURFACE 
C IN MESHING WITH A MISALIGNED REGULAR HELICAL GEAR 
C 

C NOTE 
C 

C THIS PROGRAM IS WRITTEN IN FORTRAN 77. IT CAN BE COMPILED BY V 
C COMPILER IN IBM MAINFRAME OR FORTRAN COMPILER IN VAX SYSTEM. 

C 

IMPLICIT REAL*8(A-H,0~Z) 

DOUBLE PRECISION KSIN,MU,MUO,KSIC 
DIMENSION Z (99) , ANGLE (99) , ERROR (99) , ERR (99) 

COMMON /BLOCK1/ X(ll) , Y (10) ,A(10, 10) , Y1 (10) , IPVT(10) ,WORK(10) , 

+ EPSI , DELTA , NC , NE , NDIM 

COMMON /BLOCK2/ S(4,4) ,CMM(4,4) ,CMCD(4,4) , C, R,RP , RG, Al ,HD, RL.KSIC, 
+ ALPHA , CK , SK , CB , SB , CM , SM , CA , S A , CKC , SKC , NTOOL 


C 

C DEFINE PARAMETERS USED BY PROGRAMS 
C 

C (1) IN ANG LP ARE UNIT NUMBERS ASSIGNED TO THE INPUT AND OUTPUT 
C DEVICES 

IN=5 


LP=6 

C (2) NDBUG IS USE TO CONTROL THE AUXILIARY OUTPUT FOR DEBUGGING 
NDBUG= 1 

C (3) NC IS THE UPPER LIMITATION OF REPEATATION FOR SOLVING NONLINEAR 

C EQUTIONS ; 

C EPSI IS THE CLEARANCE OF FUNCTION VALUES WHEN THE FUNCTIONS 

C IS CONSIDERED AS SOLVED (ALL FUNTIONS HAVE FORMS OF F(X)=0); 

C DELTA IS THE RELATIVE DIFFERENCE FOR TAKING DERIVATIVES 

C NC, EPSI AND DELTA MAY BE CHANGED WHEN SOLUTIONS ARE DIVERGENT 

C OR LESS ACCURATE 

NC=100 
DELTA=l.D-3 
EPSI=1.D-12 

C (4) OTHER PARAMETERS (DON'T CHANGE) 

NDIM=10 

NE=5 
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DR=DATAN ( 1 . DO) /45 . DO 

C DEFINE INPUT PARAMTERS OF PROBLEM (USE INCH AS UNIT OF LENGTH) 

C (1) PINION AND GEAR: PN=DIAMETRAL PITCH; NP=PINION TOOTH NUMBER; 

C RMPG=TOOTH NUMBER RATIO (GEAR TOOTH NO./NP); 

C KSIN=PRESSURE ANGLE IN NORMAL SECTION; 

C BETAP=HELIX ANGLE OF PINION AND GEAR; 

C HD=HEIGHT OF DEDENDUM OF PINION; 

C COE=COEFF. OF CENTRAL DISTANCE (USUALLY COE=l.) 

PN= 1 0 . DO 
NP=20 
RMPG=2 . DO 
KSIN=20.D0*DR 
BETAP= 1 5 . D0*DR 
HD=1 .DO/PN 
COE=l . OOODO 

(2) TOOL: ALPHA=HALF OF CONE VERTEX ANGLE (DEGREE) ; 

RC=RADIUS OF BOTTOM CIRCLE OF CONE 
R=RADIUS OF ARC 

MU=TILT ANGEL OF MOUNTING TOOL 

NTOOL=TOOL ID NO. (l=CONE SURFACE, 2=REVOLUTE SURFACE) 

NTOOL=2 

ALPHA=20 . DO*DR 
RC=10. 3527D0 
R=3 . OD1 

MUO=DATAN (DSIN (KSIN) *DTAN (BETAP) ) 

MU= 0 . *MUO 

(3) MISALIGNMENT: NMIS=ID NO. (l=CROSSING AXES, 2= INTERSECTING AXES); 
NG=NO. OF MISALIGNED ANGLES TO BE SIMULATED (FROM -(NG~l)/2 TO 

(NG-l)/2 TIMES GAMMAI WITH ODD NG) ; 

GAMMAI= INCREMENT OF MISALIGNED ANGLE (MINUTE) ; 

NMIS=2 
NG=2 

GAMMAI =5. DO 

C (4) OUTPUT: FEE 1= INCREMENT OF ROTATION ANGLE OF PINION (DEGREEE) 
FEEI=1.0D0*DR 
C 

C DESCRIPTION OF OUTPUT PARAMERTERS 
C 

C FEEl=ROTATION ANGLE OF PINION 

C FEE2=ROTATION ANGLE OF GEAR 

C RP=RADIUS OF PINION CONTACT POINT 

C RG=RADIUS OF GEAR CONTACT POINT 

C 

C FIND AUXILIARY VALUES FOR CALCULATION 
C 

c 

SK=DSIN (KSIN) 

CK=DC0S (KSIN) 

SB=DSIN (BETAP) 

CB=DC0S (BETAP) 

SM=DSIN (MU) 

CM=DC0S (MU) 

SA=DS IN (ALPHA) 

CA=DCOS (ALPHA) 


114 


non 


KSIC=DATAN (SK/CK/CB) 

CKC=DCOS(KSIC) 

SKC=DSIN (KSIC) 

PT=PN*CB 
RP=NP/2./PT 
RB=RP*DCOS (KSIC) 

RA=RP+1./PN 
NG=NP*RMPG 
RG=RP*RMPG 
RL=RC*CA/SA 
C CH=HD/CK/CM 

A1=RL/CA-HD/CK/CM 
C= (RP+RG) *COE 

NCOEF=360.D0*DR/FEEI/FLOAT (Nl)+0.3 

N=NC0EF*2+1 

CALL INTMAT (CMCD, 4 , 4) 

AA=RL*SA*SA/CA-HD/CK/CM 

CMCD (1 , 1) =CA*CK*CB+SA*SK*CM*CB+SA*SM*SB 

CMCD (1 , 2) =SA*CK*CB-CA*SK*CM*CB-CA*SM*SB 

CMCD (1,3) =SK*SM*CB~CM*SB 

CMCD ( 1 , 4) =-RL*SA*CK*CB-AA* (SK*CM*CB+SM*SB) 

CMCD (2, 1) =CA*SK-SA*CK*CM 

CMCD (2 , 2) =SA*SK+CA*CK*CM 

CMCD (2 , 3) =~CK*SM 

CMCD (2 , 4) =-RL*SA*SK+AA*CK*CM 

CMCD (3,1) =CA*CK*SB+SA*SK*CM*SB-SA*SM*CB 

CMCD (3 , 2) =SA*CK*SB-CA*SK*CM*SB+CA*SM*CB 

CMCD (3 , 3) =SK*SM*SB+CM*CB 

CMCD (3 , 4) =-RL*SA*CK*SB+AA* (-SK*CM*SB+SM*CB) 
CALL INTMAT (S, 4, 4) 

NGG= (NG-l)/2 
DO 505 LL=1 ,NG 

GAMMA=GAMMAI*FLOAT (LL-NGG) /60 . D0*DR 
CG=DCOS (GAMMA/60 . *DR) 

SG=DS IN (GAMMA/60 . *DR) 

IF (NMIS.EQ. 1) THEN 
WRITE (LP,500) COE, GAMMA 

500 FORMAT (1H1,///,1X, 'C =I ,F4.2, ' * (RP+RG) 

+ F5 . 1 , ' (M) ' ) 

S (1 , 1) =CG 
S(1,3)=-SG 
S (3 , 1) =SG 
S (3 , 3) =CG 
ELSE 

WRITE (LP , 501) COE, GAMMA 

501 FORMAT (1H1,///, IX, 'C=' ,F4.2, '* (RP+RG) 

+ F5 . 1 , ' (M) 1 ) 

S(2,2)=CG 
S(2,3)=-SG 
S(3,2)*SG 
S(3,3)-CG 
END IF 

DO 205 L-1,2 
DO 5 1=1, NE 


CROSSING ANGLE= ' , 


INTERSECTING ANGLE= ' 
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5 X(I)=0.D0 

IF (NTOOL.EQ. 1) X(5)=RL/CA~HD/CK/CM 
C WRITE (6,1100) (X( JK) , JK=1 , 5) ,RL, CA, HD, CK, CM 

Cl 100 FORMAT (IX, '#####' ,5F15.7/7X,5F15. 7) 

DO 15 1=1, N 

X (7) =FEEI*FLOAT (I~ (N+l) /2) 

IF (L.EQ.l) X(7)=0.D0 
C X(5)=DARSIN(RP*X(7)/SK/R) 

X (2) =X (7) /RMPG 
CALL NONLIN 
X(8)=X(2)+X(3) 

IF (L.EQ.l) THEN 
XIN=X(8) 

WRITE (LP , 10) 

10 FORMAT (////8X, 'FEEl(D) ' ,8X, 'FEE2(D) ' ,8X, 'K-ERROR(S) ' ,5X, 

+ ' RP ’ , 13X, ' RG’ , F15 . 7/) 

GO TO 205 
END IF 

X(8) =X (8) -XIN 
X(7)=X(7)/DR 
X(8)=X(8)/DR 

X (9) = (X (8) -X ( 7) /RMPG) *3600 . DO 
Z(I)=X(8) 

ERR (I) =X (9) 

WRITE (LP , 20) (X(J),J=7,11) 

20 FORMAT (1X.5F15 . 7 , F15 . 7) 

WRITE (LP , 30) (X(J) , J=1 , 6) 

30 FORMAT (IX, '#####' ,6F15.7) 

WRITE (LP,30) XP , YP ,ZP ,XG, YG 
15 CONTINUE 
NT=N-NCOEF 

FII=FEEI/DR/2 . DO*FLOAT (NCOEF) 

WRITE (LP , 80) 

80 FORMAT (//, ' FIND THE WORKING RANGE FOR ONE TOOTH: ' , Fl5 . 7/) 
DO 55 1=1, NT 

X (7) =FEEI*FLOAT (l~ (N+l) /2) /DR/2 . DO 
X(8)=X(7)+FII 
KK=I+NCOEF 
ANGLE (I) =Z (KK) -Z (I) 

ERROR (I) = (ANGLE (I) -FIl) *3600. DO 
WRITE (LP ,60) X (7) ,X (8), ANGLE (I), ERROR (I) 

60 FORMAT (IX, ' ( ' , F7. 2 , ' ' , F7 . 2, ’ ) : ' , F15 . 7 , F15 . 7) 

55 CONTINUE 
DO 95 1=1, NT 
ATEMP2=ERR0R(I) 

IF (I.NE.1) THEN 

IF (ATEMP1*ATEMP2.LE.0.D0) GOTO 105 
END IF 

ATEMP1=ATEMP2 
95 CONTINUE 

WRITE (LP , 160) 

160 FORMAT (//IX, 'MESHING IS DISCONTINUOUS') 

105 IF (DABS (ATEMP1) .LT. DABS (ATEMP2)) 1=1-1 
EMAX=0. 
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EMIN=0. 

DO 135 J=1 ,NCOEF 
KS=I+J-1 
ET=ERR (KS) 

IF (ET . LT. EMIN) EMIN=ET 
IF (ET.GT.EMAX) EMAX=ET 
135 CONTINUE 

ET=EMAX-EMIN 

KK=I+NCOEF 

WRITE (LP, 170) Z(l) ,Z(KK) ,ET 

170 FORMAT (//IX, 'WORKING RANGE FOR ONE TOOTH: ',F7.2,' ',F7.2/ 

+ IX, 'THE MAXIMUM KINEMATIC ERROR: '.F15.7,' (S) ' , 12) 

205 CONTINUE 
505 CONTINUE 
STOP 
END 


C 

C 


SUBROUTINE FUNC 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,MUO,KSIC 

COMMON /BLOCK1/ X (11) , Y (10) , A(10, 10) , Yl (10) , IPVT (10) , WORK(IO) , 
h EPSI .DELTA, NC,NE,NDIM 

COMMON /BLOCK2/ S(4,4) ,CMM(4,4) ,CMCD(4,4) ,C,R,RP ,RG,Al ,HD,RL,KSIC, 
- ALPHA , CK , SK , CB , SB , CM , SM , CA , S A , CKC , SKC , NTOOL 

DIMENSION RD(4) ,RND(4) ,RC(4) ,RNC(4) ,RF(4) 

CFEE=DCOS (X (3) ) 

SFEE=DSIN (X (3) ) 

CT“DCOS (X (4) ) 

ST=DSIN (X (4) ) 

CF-DCOS (X (7) ) 

SF=DSIN (X(7) ) 

IF (NTOOL. EQ.l) THEN 

RD(1)=X(5)*CT*SA 

RD(2)=RL-X(5)*CA 

RD (3) =X (5) *ST*SA 

RND(l) =CT*CA 

RND(2)=SA 

RND(3)=ST*CA 

ELSE 

CAL=DCOS (ALPHA+X (5) ) 

SAL=DSIN (ALPHA+X (5) ) 

SL2=DSIN(X(5)/2.) 

CAL2=DCOS (ALPHA+X (5) /2 . ) 

SAL2=DS IN- (ALPHA+X (5) /2 . ) 

RD(1)= (Al*SA-2. *R*SL2*SAL2) *CT 

RD (2) =HD/CK/CM*CA+2 . *R*SL2*CAL2 

RD (3) = (Al*SA-2 . *R*SL2*SAL2) *ST 

RND(l) =CAL*CT 

RND(2)=SAL 

RND (3) =CAL*ST 

END IF 

RD(4) =1 . 

CALL MATMUL (CMCD , RD , RC , 4 , 4) 


117 


non nn nnnn 


CALL MATMUL(CMCD,RND,RNC,3,3) 

X (6) = (RC (1) -RC (2) *RNC (1) /RNC (2) ) /RP 
RF (1)=RC (1) -RP*X (6) 

RF (2) =RC (2) +RP 
RF (3) =RC (3) 

RF (4) =1 . 

CFPF1=DC0S (X (6) +X (7) ) 

SFPFl=DSIN (X (6) + X (7) ) 

XPF=RF (1) *CFPFl+RF (2) *SFPFl 
YPF=RF (2) *CFPF1-RF (1) *SFPF1 
ZPF=RF (3) 

XNPF=RNC (1) *CFPF1+RNC (2) *SFPF1 
YNPF=RNC (2) *CFPFl-RNC (1) *SFPF1 
ZNPF=RNC (3) 

CF2FG=DC0S (X (3) ) 

SF2FG=DSIN (X (3) ) 

A2=-X (1) *SB+RG*X (2) 

RFG1=CK*CK/CKC*DC0S (X (3) +KSIC) *A2+RG*SF2FG 

RFG2=CK*CK/CKC*DSIN (X (3) +KS IC) *A2~RG*CF2FG 

RFG3=X (1) * (CB+SK*SK*SB*SB/CB) -RG*X (2) *SK*SK*SB/CB 

RNFG1=CK*CB*CF2FG-SK*SF2FG 

RNFG2=CK*CB*SF2FG+SK*CF2FG 

RNFG3=CK*SB 

XGF=RFG1*S (1 , 1) +RFG2*S (1 , 2) +RFG3*S (1,3) 

YGF=RFG1*S (2 , 1) +RFG2*S (2 , 2) +RFG3*S (2 , 3) 

ZGF=RFG1*S (3 , 1) +RFG2*S (3 , 2) +RFG3*S (3 , 3) 

XNGF=RNFG1*S (1 , 1) +RNFG2*S (1 , 2) +RNFG3*S (1 , 3) 

YNGF=RNFG1*S (2,1) +RNFG2*S (2 , 2) +RNFG3*S (2 , 3) 

ZNGF=RNFG1*S (3,1) +RNFG2*S (3 , 2) +RNFG3*S (3 , 3) 

WRITE (6,100) X(l) ,X(2) ,X(3) ,X(4) ,X(5) 

WRITE (6,100) XPF , YPF , ZPF , XGF , YGF , ZGF 
100 FORMAT (IX, 'ZIVL'L' ,8F15.7) 

WRITE (6,100) XNP F , YNP F , ZNP F , XNGF , YNGF , ZNGF 

Y(1)=XPF-XGF 

Y(2)=YPF-YGF-C 

Y(3)=ZPF-ZGF 

Y(4)=YNPF-YNGF 

Y(5)=ZNPF-ZNGF 

X (10) -DSQRT (XPF*XPF+YPF*YPF) 

X ( 1 1) =DSQRT (XGF*XGF+YGF*YGF) 

WRITE (6,20) (Y(II) ,11=1,5) 

20 FORMAT (IX, '$$$$', 6F15 . 7) 

RETURN 

END 

SUBROUTINE INTMAT (A,N,M) 

THIS SUBROUTINE IS USED TO INITIATE THE MATRIX, WITH UNIT DIAGONAL 
ELEMENTS AND NULL OTHER ELEMENTS 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A (4, 4) 

DO 5 1=1, N 
DO 5 J=1,M 
A (I , J) =0. 

IF (I.EQ.J) A(I , J) =1 . 
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5 CONTINUE 
RETURN 
END 

SUBROUTINE MATMUL (CMM, A, B , N,M) 

THIS SUBROUTINE IS USED TO MULTIPLY THE MATRIX CMCD(N*M) BY THE MATRIX 
A(M*1). THE RESULT IS STORED IN THE MATRIX B(N*1) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION CMM (A , 4) , A (4) , B (4) 

DO 5 1=1, N 
5 B (l)=0. 

DO 15 1=1, N 
DO 15 J=1,M 

15 B(I)=B(I)+CMM(I,J)*A(J) 

RETURN 

END 


SUBROUTINE NONLIN 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON /BLOCK1/ X (1 1) , Y (10) , A (10 , 10) , Y1 (10) , IPVT (10) , WORK (10) , 
+ EPSI, DELTA, NC,NE,NDIM 

DO 5 1=1, NC 
CALL FUNC 

WRITE (6,10) I, (X(J) ,Y(J) , J-1,5) 

10 FORMAT (IX, '***' ,I5/5(1X,2D15.7/)) 

DO 15 J=1 ,NE 

IF (DABS (Y (J) ) .GT. EPSI) GO TO 25 
15 CONTINUE 
GO TO 105 
25 DO 35 J=1 ,NE 
35 Y1(J)=Y(J) 

DO 45 J=1 , NE 
DIFF=DABS (X ( J) ) *DELTA 
IF (X(J) .EQ.O.DO) DIFF=DELTA 
XMAM=X(J) 

X ( J) =X ( J) -DIFF 
CALL FUNC 
X(J) =XMAM 
DO 55 K=1 ,NE 

55 A(K, J) = (Y1(K)-Y(K))/DIFF 
45 CONTINUE 
DO 65 J=1 ,NE 
65 Y(J)=-Y1(J) 

CALL DECOMP (NDIM,NE,A,COND, IPVT, WORK) 

CALL SOLVE (NDIM,NE,A, Y, IPVT) 

DO 75 J=1,NE 
75 X(J)=X(J)+Y(J) 

5 CONTINUE 
105 RETURN 
END 
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SUBROUTINE DECOMP (NDIM,N,A,COND, IPVT,WORK) 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A(NDIM,N) ,WORK(N) ,IPVT(N) 

C 

C DECOMPOSES AREAL MATRIX BY GAUSSIAN ELIMINATION, 

C AND ESTIMATES THE CONDITION OF THE MATRIX. 

C 

C -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS", BY G. E. FORSYTHE, 
C M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

C 

C USE SUBROUTINE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEM. 

C 

C INPUT. . 

C 

C NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A 

C N ORDER OF THE MATRIX 

C A MATRIX TO BE TRIANGULAR I ZED 

C 

C OUTPUT . . 

C 

C A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PREMUTED 
C VERSION OF A LOWER TRIANGULAR MATRIX I~L SO THAT 

C (PERMUTATION MATRIX) *A=L*U 

C 

C COND = AN ESTIMATE OF THE CONDITION OF A. 

C FOR THE LINEAR SYSTEM A*X = B , CHANGES IN A AND B 

C MAY CAUSE CHANGES COND TIMES AS LARGE IN X. 

C IF COND+l.O .EQ. COND , A IS SINGULAR TO WORKING 

C PRECISION. COND IS SET TO 1 . OD+32 IF EXACT 

C SINGULARITY IS DETECTED. 

C 

C IPVT = THE PIVOT VECTOR 

C IPVT (K) = THE INDEX OF THE K-TH PIVOT ROW 

C IPVT (N) = (-1)** (NUMBER OF INTERCHANGES) 

C 

C WORK SPACE.. THE VECTOR WORK MUST BE DECLARED AND INCLUDED 
C IN THE CALL. ITS INPUT CONTENTS ARE IGNORED. 

C ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. 

C 

C THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY 
C DET(A) = IPVT (N) * A(1 , 1) * A(2,2) * ... * A(N,N) . 

C 

IPVT (N) =1 

IF (N.EQ.l) GO TO 150 
NM1=N-1 

C COMPUTE THE- 1-NORM OF A . 

ANORM=O.DO 
DO 20 J=1,N 
T=0.D0 
DO 10 1=1, N 
10 T=T+DABS(A(I, J)) 
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IF (T.GT.ANORM) ANORM=T 
20 CONTINUE 

DO GAUSSIAN ELIMINATION WITH PARTIAL 
PIVOTING. 

DO 70 K=l,NMl 
KP1=K+1 

FIND THE PIVOT. 

M=K 

DO 30 I=KP1,N 

IF (DABS (A(l ,K) ) .GT.DABS (A(M,K) ) ) M=I 
30 CONTINUE 
IPVT(K) =M 

IF (M.NE.K) IPVT(N)— IPVT(N) 

T=A (M,K) 

A (M , K) =A(K,K) 

a(k,k)=t 

SKIP THE ELIMINATION STEP IF PIVOT IS ZERO. 
IF (T.EQ.O.DO) GO TO 70 

COMPUTE THE MULTIPLIERS. 

DO 40 I=KP1,N 
40 A(I,K)=-A(I,K)/T 

C INTERCHANGE AND ELIMINATE BY COLUMNS. 

DO 60 J=KPl,N 
T=A (M, J) 

a(m, j) =a(k, j) 

A(K, J)=T 

IF (T.EQ.O.DO) GO TO 60 
DO 50 I=KPl,N 

50 A(I, J)=A(I, J)+A(I,K)*T 
60 CONTINUE 
70 CONTINUE 
C 

C COND = (1-NORM OF A)* (AN ESTIMATE OF THE 1-NORM OF A-INVERSE) 

C THE ESTIMATE IS OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE 
C SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS 
C OF EQUATIONS, (A- TRANSPOSE) *Y = E AND A*Z = Y WHERE E 
C IS A VECTOR OF +1 OR -1 COMPONENTS CHOSEN TO CAUSS GROWTH IN Y. 

C ESTIMATE = (1-NORM OF Z) / (1-NORM OF Y) 

C 

C SOLVE (A- TRANSPOSE) *Y = E . 

DO 100 K=1,N 
T=0.D0 

IF (K.EQ.l) GO TO 90 
KM1=K-1 
DO 80 1=1, KM1 
80 T=T+A(I,K)*W0RK(I) 

90 EK= 1 . DO 

IF (T.LT.O.DO) EK=-1.D0 
IF (A(K,K) .EQ.O.DO) GO TO 160 
100 WORK (K) =- (EK+T) /A (K , K) 

DO 120 KB=1 ,NMl 
K=N-KB 
T=0.D0 
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KP1=K+1 

DO 110 I=KPl , N 
110 T=T+A(l,K)*WORK(K) 

WORK (K) =T 
M=IPVT(K) 

IF (M.EQ.K) GO TO 120 
T=WORK (M) 

WORK (M) =WORK (K) 

WORK (K) =T 
120 CONTINUE 
C 

YNORM=O.DO 
DO 130 1=1, N 

130 YNORM=YNORM+DABS (WORK (I) ) 

C 

C SOLVE A*Z = Y 

CALL SOLVE (NDIM,N,A, WORK, IPVT) 

C 

ZNORM=0 . DO 
DO 140 1=1, N 

140 ZNORM=ZNORM+DABS (WORK (i) ) 

C 

C ESTIMATE THE CONDITION. 

COND= ANORM*ZNORM/ YNORM 

IF (COND.LT. l.DO) COND=l.DO 

RETURN 

C 1-BY-l CASE.. 

150 COND= 1 . DO 

IF ( A ( 1 , 1 ) . NE . 0 . DO) RETURN 
C 

C EXACT SINGULARITY 

160 C0ND=1.0D32 
RETURN 
END 

SUBROUTINE SOLVE (NDIM, N , A, B , IPVT) 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A (NDIM, N) , B (N) , IPVT (N) 

C 

C SOLVES A LINEAR SYSTEM, A*X = B 

C DO NOT SOLVE THE SYSTEM IF DECOMP HAS DETECTED SINGULARITY. 

C 

C -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS-, BY G. E. FORSYTHE, 
C M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

C 

C INPUT.. 

C 

C NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A 
C N ORDER OF MATRIX 

C A TRIANGULAR I ZED MATRIX OBTAINED FROM SUBROUTINE DECOMP 

C B RIGHT HAND SIDE VECTOR 

C IPVT = PIVOT VECTOR OBTAINED FROM DECOMP 

C 

C OUTPUT.. 
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B = SOLUTION VECTOR, X 


DO THE FORWARD ELIMINATION. 

IF (N.EQ.l) GO TO 50 
NM1-N-1 
DO 20 K=1 ,NMl 
KP1=K+1 
M=IPVT(K) 

T=B (M) 

B (M) =B (K) 

B (K) =T 

DO 10 I-KPl.N 
10 B(I)=B(I)+A(I,K)*T 

20 CONTINUE 

C NOW DO THE BACK SUBSTITUTION. 

DO 40 KB=1 , NMl 
KMl=N-KB 
K=KM1+1 

B(K)=B(K)/A(K,K) 

T— B (K) 

DO 30 1=1, KM1 
30 B(l)=B(l)+A(I,K)*T 

40 CONTINUE 
50 B (1) =B (1) /A (1 , 1) 

RETURN 

END 

C 

C 
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******************************************************* 
* * 

* * 


* PROGRAM V * 

* TRANSMISSION ERRORS OF CROWNED HELICAL PINION IN * 

* MESHING WITH REGULAR HELICAL GEAR WITH PREDESIGNED * 

* TRANSMISSION ERRORS AND CONTACT PATH * 


* * 

* AUTHORS: FAYDOR LITVIN * 

* JIAO ZHANG * 

* * 


* * 
******************************************************* 


C 

C PURPOSE 

C 

C 

C THIS PROGRAM IS USED TO CALCULATE THE TRANSMISSION ERRORS OF A 
C CROWNED HELICAL PINION AND A HELICAL GEAR WHEN THEIR AXES ARE 
C MOUNTED WITH SOME ERRORS. THE MAGNITUDE OF TRANSMISSION ERRORS 
C AND CONTACT PATH ARE PRE-DESIGNED. 

C 

C NOTE 
C 

C THIS PROGRAM IS WRITTEN IN FORTRAN 77. IT CAN BE COMPILED BY V 
C COMPILER IN IBM MAINFRAME OR FORTRAN COMPILER IN VAX SYSTEM. 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,KK,KSIC 

DIMENSION Z (99) .ANGLE (99) .ERROR (99) .ERR (99) ,W(99) , RPR (99) ,RGR(99) , 
+ ZP (99) ,ZG(99) 

COMMON /BLOCK1/ X (13) , Y (10) , A (10, 10) , Y1 (10) , IPVT (10) , WORK (10) , 

+ EPSI, DELTA, NC.NE.NDIM 

COMMON /BLOCK2/ S(4,4) ,AT(5) .C.RP.RG.CK.SK.CB.SB.CKC.KSIC.TB, 

+ COEG 1 , COEG2 , RMGP , AA , DR , KK , R 

C 

C DEFINE PARAMETERS USED BY PROGRAMS 
C 

C (1) IN ANG LP ARE UNIT NUMBERS ASSIGNED TO THE INPUT AND OUTPUT 
C DEVICES 

IN=5 
LP=6 

C (2) NDBUG IS USE TO CONTROL THE AUXILIARY OUTPUT FOR DEBUGGING 
NDBUG= 1 

C (3) NC IS THE UPPER LIMITATION OF REPEATATION FOR SOLVING NONLINEAR 
C EQUTIONS; 

C EPSI IS THE CLEARANCE OF FUNCTION VALUES WHEN THE FUNCTIONS 

C IS CONSIDERED AS SOLVED (ALL FUNTIONS HAVE FORMS OF F(X)=0); 

C DELTA IS THE RELATIVE DIFFERENCE FOR TAKING DERIVATIVES 

C NC, EPSI AND DELTA MAY BE CHANGED WHEN SOLUTIONS ARE DIVERGENT 

C OR LESS ACCURATE 

NC=100 
DELTA=l.D-3 
EPSI=1.D-10 
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C (4) OTHER PARAMETERS (DON'T CHANGE) 

NDIM=10 

NE=5 

DR=DATAN ( 1 . DO) /45 . DO 

C DEFINE INPUT PARAMTERS OF PROBLEM (USE INCH AS UNIT OF LENGTH) 

C (1) PINION AND GEAR: PN=DIAMETRAL PITCH; NP=PINION TOOTH NUMBER; 

C RMPG=TOOTH NUMBER RATIO (GEAR TOOTH NO./NP); 

C KSIN=PRESSURE ANGLE IN NORMAL SECTION; 

C BETAP=HELIX ANGLE OF PINION AND GEAR; 

C COE=COEFF . OF CENTRAL DISTANCE (USUALLY COE=l.) 

PN=2.D0 

NP=12 

RMPG=94./12. 

KSIN=30.D0*DR 
BETAP=15 . DO^DR 
COE=l . OOODO 

C (2) PREDESIGN TRANSMISSION ERRORS AND CONTACT PATH: 

C AA=LEVEL OF PREDESIGNED PARABOLIC TRNASMISSION ERRORS IN SECOND 

C KK=COEFFICIENT OF CONTACT PATH DIRECTION (l.D~3=CROSS TOOTH 

C SURFACE; 5 . D6=ALONG TOOTH SURFACE) 

C R=RADIUS OF ARC ATTACHED TO CHOSEN CONTACT PATH 

AA=25.D0 
KK=5 . D6 
R=0 . 3584D0 

(3) MISALIGNMENT: NMIS=ID NO. (l=CROSSING AXES, 2= INTERSECTING AXES); 
NG=NO. OF MISALIGNED ANGLES TO BE SIMULATED (FROM -(NG-l)/2 TO 
(NG-l)/2 TIMES GAMMAI WITH ODD NG) ; 

GAMMAI= INCREMENT OF MISALIGNED ANGLE (MINUTE) ; 

NMIS=1 
NG=3 

GAMMAI =3. DO 

C (4) OUTPUT: FEEI=INCREMENT OF ROTATION ANGLE OF PINION (DEGREEE) 
FEEI=1.0D0*DR 
C 

C DESCRIPTION OF OUTPUT PARAMERTERS 
C 

C FEEl'ROTATION ANGLE OF PINION 

C FEE2=ROTATION ANGLE OF GEAR 

C RP=RADIUS OF PINION CONTACT POINT 

C RG=RADIUS OF GEAR CONTACT POINT 

C ZP=LENGTH OF PINION CONTACT POINT FROM MIDDLE SECTION 

C ZG=LENGTH OF GEAR CONTACT POINT FROM MIDDLE SECTION 
C 

C FIND AUXILIARY VALUES FOR CALCULATION 
C 

C DEFINE USEFUL CONSTANTS AND PARAMETERS FOR PINION AND GEAR 
RMGP=1./RMPG 
NG=NP*RMPG+0 . 5 

AA=AA* (2/3600 . *DR* (NP/DR/180 . ) **2 ) 

SK=DSIN (KSIN) 

CK=DC0S (KSIN) 

SB=DSIN (BETA) 

CB=DCOS (BETA) 

TB=SB/CB 
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KSIC'DATAN (SK/CK/CB) 

CKC=DCOS (KSIC) 

SKC=DSIN (KSIC) 

PT=PN*CB 

RP=NP/2./PT 

RPB=RP*CKC 

RPA=RP+1 . O/PN 

RG=NG/2./PT 

RGB=RG*CKC 

RGA=RG+ 1 . O/PN 

WRITE (LP,56) RP , RPB , RPA , RG , RGB , RGA , AA 
56 FORMAT (IX, '&&&&&&', ' RP= ' , F15 . 7, 5X, 'RPB= 1 ,F15.7,5X, ' RPA= ' , F15 . 7/ 
+ IX, * &&&&&& ' , 'RG= ' , F15 . 7, 5X, ’RGB*' ,F15.7,5X, ' RGA= ' ,F15.7) 

CONl= (1 . +SK*SK*TB*TB) 

AT (1) = (CK**4/CKC**2*KK**2~SK**4*TB**2) *RG*RG 

AT (2) = (2 . *SK*SK*SB*CONl-2 . *KK**2*SB*CK**4/CKC**2) *RG 

AT (3) = (-2 . *CK*CK/CKC*SKC*KK**2-2 . *KK*SK*SK*TB) *RG*RG 

AT (4) = (2 . *SB*CK*CK/CKC*SKC*KK**2+2 . *KK*CB’''CONl) *RG 

AT (5) =CK**4/CKC**2*SB**2*KK**2-CB**2*CONl*CONl 

C= (RP+RG) *COE 

CALL INTMAT(S,4,4) 

NGG= (NG-l)/2 

DO 505 LL=1 ,NG 

GAMMA= GAMMA I * F LO AT (LL-NGG) 

CG=DCOS (GAMMA/60. *DR) 

SG=DS IN (GAMMA/60. *DR) 

IF (NMIS.EQ. 1) THEN 
WRITE (LP , 500) COE, GAMMA 

500 FORMAT (1H1,///,1X, 1 C— ' ,F9.4, ' * (RP+RG) CROSSING ANGLE= ' , 

+ F8 . 4 , ' (M) ' ) 

S (1 , 1) =CG 
S(1,3)=-SG 
S (3 , 1) =SG 
S(3,3)=CG 
ELSE 

WRITE (LP,501) COE, GAMMA 

501 FORMAT (1H1,///, IX, ' C- ' ,F4.2, '* (RP+RG) INTERSECTING ANGLE= ' , 

+ F5 . 1 , ' (M) ' ) 

S (2, 2) =CG 
S(2,3)=-SG 
S (3 , 2) =SG 
S (3 , 3) =CG 
END IF 

NCOEF=IDINT (360 . *DR/ FEE I /FLOAT (NP) +0 . 5) 

N=2*NCOEF 
NHALF= (N+l)/2 
XIN=0. 

DO 205 L=1 , 2 
LSGN= (-1)**L 
DO 15 1=1 , NHALF 
LI=NHALF+LSGN*(I-1) 

IF (L.EQ.l) NMIN=LI 
IF (L.EQ.2) NMAX=LI 
X (7) =LSGN*FEEI*FLOAT (i-l) 
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X (4) =X(7) 

X (2) =X (7) /RMPG 
X(3)=0. 

X (5) = (- (AT (2) *X (2) +AT (A) ) +DSQRT ( (AT (2) *X (2) +AT (4) ) **2~4 . *AT (5) 
+ * (AT (3) *X (2) +AT (1) *X (2) *X (2) ) ) ) /2 . /AT (5) 

X (1) =0. 

CALL NONLIN 
X(8)=X(1)+X(2) 

C FIND INITIAL VALUE OF X(8) 

IF (L.EQ.l.AND.I.EQ.l) THEN 
XIN=X(8) 

WRITE (LP,51) XIN 
51 FORMAT (IX, 'XIN=' ,D15.7) 

GO TO 15 
END IF 
C 

X (8) =X (8) -XIN 
W(LI)=X(7)/DR 
Z(LI)=X(8)/DR 

ERR (LI) = (X (8) -X (7) /RMPG) *3600 . DO/DR 
RPR (Li) =X (10) 

RGR(LI)=X(11) 

ZP (Li) =X (12) 

ZG (LI) =X (13) 

C WRITE (LP , 20) W(LI) ,Z(LI) , ERR (LI) , RPR (LI) ,RGR(LI) 

15 CONTINUE 
205 CONTINUE 

WRITE (LP, 10) 

10 FORMAT (////8X, 'FEEl(D) ' ,8X, ’FEE2(D) ' ,8X, 'K-ERROR(S) ' ,5X, 

+ ' RP ' , 13X, ' ZP ' , 13X, 'RG' , 13X, ' ZG ' , F15 . 7/) 

DO 25 I=NMIN,NMAX 

WRITE (LP , 20) W(l) ,Z(l) , ERR (i) , RPR (I) , ZP (I) , RGR (I) , ZG (i) 

20 FORMAT (IX, 2F15 . 7 , F12. A , 4F15 . 7) 

25 CONTINUE 

NT=NMAX-NCOEF 

FII=FEEl/DR/RMPG*FLOAT (NCOEF) 

WRITE (LP,80) 

80 FORMAT (//,' FIND THE WORKING RANGE FOR ONE TOOTH: ' , F15 . 7/) 

DO 55 I=NMIN,NT 

X (7) =FEEI*FLOAT (l~ (N+l) /2) /RMPG/DR 
X(8)=X(7)+FII 
KK=I+NCOEF 
ANGLE (I) =Z (KK) ~Z (I) 

ERROR (I) = (ANGLE (I) -FII) *3600 . DO 
WRITE (LP , 60) X(7) ,X(8) , ANGLE ( I) , ERROR (I) 

60 FORMAT (IX, ' ( ' , F7. 2 , ' ' , F7.2, ' ) : ' , F15. 7, F15. 7) 

55 CONTINUE 

DO 95 I=NMIN,NT 
ATEMP2=ERR0R(I) 

IF (I.NE.NMIN) THEN 
IF (ATEMP 1 *ATEMP2 . LE . 0 . DO) GOTO 105 
END IF 

ATEMP 1= ATEMP 2 
95 CONTINUE 
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WRITE (LP, 160) 

160 FORMAT (//IX, 'MESHING IS DISCONTINUOUS') 

GO TO 505 

105 IF (DABS (ATEMP1) .LT. DABS (ATEMP2)) 1=1-1 
EMAX=0. 

EMIN=0 . 

NTEMP-NCOEF+1 
DO 135 J=1 , NTEMP 
KS-I+J-1 
ET=ERR (KS) 

IF (ET.LT.EMIN) EMIN=ET 
IF (ET.GT.EMAX) EMAX=ET 
135 CONTINUE 

ET=EMAX-EMIN 

KK-I+NCOEF 

WRITE (LP , 170) Z (I) , Z (KK) , ET 

170 FORMAT (//IX, 'WORKING RANGE FOR ONE TOOTH: \F7.2,' ',¥ 1 . 2 / 

+ IX, 'THE MAXIMUM KINEMATIC ERROR: '.F12.A,' (S) ' , 12) 

505 CONTINUE 
STOP 
END 

SUBROUTINE FUNC 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,KK,KSIC 

COMMON /BLOCK1/ X (13) , Y (10) , A (10 , 10) , Y1 (10) , IPVT (10) , WORK(IO) , 
+ EPSI ,DELTA,NC,NE,NDIM 

COMMON /BLOCK2/ S(A,A) ,AT(5) ,C,RP,RG,CK,SK,CB,SB,CKC,KSIC,TB, 

+ COEG1 ,COEG2,RMGP,AA,DR,KK,R 

*********>V**********5Wc**************************** 

20 TEMP-CKC* ( 1 . -AA*X (A) / (RMGP+ 1 . ) ) 

IF (DABS (TEMP) . GT . 1 . DO) THEN 
X (A) =X (7) 

GOTO 20 
END IF 

WRITE (6,220) TEMP , X (A) , AA 

220 FORMAT (IX, 'TEMP-' ,E15.7, 'X(A)=' ,El5.7, 'AA-' ,F15.7) 
RLAM-KSIC-DARCOS (TEMP) 

DLAMDF—CKC*AA/ (RMGP+1 . ) /DSQRT (1 .~TEMP*TEMP) 

FEEGP-X (A) *RMGP-AA/2 . *X (A) *X (A) +RLAM 
DFGDFP=RMGP-AA*X (A) +DLAMDF 
CONA- AT (5) 

CONB-AT (2) *FEEGP+AT (A) 

CONC-AT (3) *FEEGP+AT (1) *FEEGP*FEEGP 
COND=CONB*CONB-A . *CONA*CONC 
IF (COND.LT.O.) THEN 
X (A) -X (7) 

GOTO 20 
END IF 

COND-DSQRT (COND) 

TP= (-CONB+COND) /2 . /CONA 
IF (COND. EQ. 0.0) THEN 
DTPDFG— AT (2) /2 . /CONA 
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ELSE 

DTPDFG= (-AT (2) - (AT (2) *C0NB-2 . *CONA* (AT (3) +2 . *AT ( 1) *FEEGP) ) /COND) 
+ /2./CONA 


END IF 


C *************************************************** 

C X(4)=FEEP; X (7) =FEE1 
APl=X (4) -X(7) 

AP2=AP 1-RLAM 
AP3=AP2+KSIC 
AP4=X (4) -RLAM+KSIC 
DAP4DF= 1 . -DLAMDF 
CAP4=DCOS (AP4) 

SAP4=DSIN (AP4) 

AP=-TP*SB+RG*FEEGP 

DAPDF= (~SB*DTPDFG+RG) *DFGDFP 

MU=X (4) +KSIC-RLAM 

DMUDFP=1 . +DLAMDF 

ALP1=R* (DSIN (MU+X (3) ) -DSIN (MU) ) 

ALP2=R* (DCOS (MU+X (3) ) -DCOS (MU) ) 

XPF=CK*CK/CKC*DCOS (AP3) *AP+RG*DSIN (AP2) - (RP+RG) *DSIN (API) +ALP2 
YPF=CK*CK/CKC*DSIN (AP3) *AP~RG*DCOS (AP2) + (RP+RG) *DCOS (API) +ALP1 
ZPF=TP*CB*CK*CK/CKC/CKC-SK*SK*TB*RG*FEEGP 

Q'k'k'k'k'k'kjz’k'k'k'k'k'ick'k'kjck'fc’k'kjt'klck'jtirkjck'jrirk'kic'k&’it’k-k’fc'is'k'kjt'jtirkJck ******** 


DXDA — DSIN (MU+X (3) ) 

DYDA=DCOS (MU+X (3) ) 

DXDF=CK*CK/CKC ,V (CAP4*DAPDF-SAP4*AP*DAP4DF) +RG*DCOS (X (4) -RLAM) 
+ * ( 1 . -DLAMDF) - (RP+RG) *DCOS (X (4) ) -ALP 1 *DMUDFP 

DYDF=CK*CK/CKC* (SAP4*DAPDF+CAP4*AP*DAP4DF) +RG*DSIN (X (4) -RLAM) 
+ * (1 . -DLAMDF) - (RP+RG) *DSIN (X (4) ) +ALP2*DMUDFP 

DZDF=CB*CK*CK/CKC/CKC*DTPDFG*DFGDFP-SK*SK*TB*RG*DFGDFP 
XNP=DYDA*DZDF 
YNP=-DXDA*DZDF 
ZNP=DXDA*DYDF-DYDA*DXDF 
RMN=DSQRT (XNP*XNP+YNP*YNP+ZNP*ZNP) 

C WRITE (6 , 130) DXDA , DZDA , DXDF , DYDF , DZDF , XNP , YNP , ZNP , RMN , TP , AP 
C 130 F0RMAT(1X, '$$$' , 5F15 . 7/4X.5F15 . 7) 

XNP=XNP/RMN 

YNP=YNP/RMN 

ZNP=ZNP/RMN 

QjrkjtftJc'kJijc'k'k'kjc'kicjtickjcit'kJt'k'kTfc'kic'k'kjc'k'k’kicic'kjc'k'k'k'k'itjtjcfc'k'icjckjt'k'k'k'kjc’k'k'Jck 


XNPF=XNP*DCOS (X (7) ) +YNP*DSIN (X (7) ) 
YNPF=YNP*DCOS (X (7) ) -XNP*DSIN (X (7) ) 
ZNPF=ZNP 
CF2FG=DCOS (X(l)) 

SF2FG=DSIN(X(1)) 

CF2FGK=DCOS (X (1) +KSIC) 

SF2FGK=DSIN (X (1) +KSIC) 

AG=-X(5)*SB+RG*X(2) 

RFG1=CK*CK/CKC*CF2FGK*AG+RG*SF2FG 

RFG2=CK*CK/CKC*SF2FGK*AG-RG*CF2FG 

RFG3=X (5) *CB-AG*SK*SK*TB 

RNFG1=CK*CB*CF2FG-SK*SF2FG 

RNFG2=CK*CB*SF2FG+SK*CF2FG 

RNFG3=CK*SB 
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XGF=RFG1*S (1 , 1) +RFG2*S (1 , 2) +RFG3*S (1 , 3) 

YGF=RFG1*S (2,1) +RFG2*S (2 , 2) +RFG3*S (2 , 3) 

ZGF=RFG1*S (3,1) +RFG2*S (3 , 2) +RFG3*S (3 , 3) 

XNGF=RNFG1*S (1 , 1) +RNFG2*S (1,2) +RNFG3*S (1,3) 

YNGF=RNFG1*S (2,1) +RNFG2*S (2 , 2) +RNFG3*S (2 , 3) 

ZNGF=RNFG1*S (3 , 1) +RNFG2*S (3 , 2) +RNFG3*S (3 , 3) 

WRITE (6,100) X(l) ,X(2) ,X(3) ,X(4) ,X(5) 

WRITE (6,100) XPF , YPF , ZPF , XGF , YGF , ZGF 
100 FORMAT (IX, 'IV/.1V ,8E15.7) 

WRITE (6,100) XNPF , YNPF , ZNPF , XNGF , YNGF , ZNGF 

Y(1)=XPF-XGF 

Y(2)=YPF-YGF-C 

Y(3)=ZPF-ZGF 

Y(5)=XNPF-XNGF 

Y (4) =ZNPF-ZNGF 

X (10) =DSQRT (XPF*XPF+YPF*YPF) 

X (1 1) =DSQRT (XGF*XGF+YGF*YGF) 

X (12) =ZPF 
X(13)=ZGF 

WRITE (6,20) (Y(II) ,11=1,5) 

20 FORMAT (IX, '$$$$’, 6F15 . 7) 

RETURN 

END 

SUBROUTINE INTMAT (A,N,M) 

THIS SUBROUTINE IS USED TO INITIATE THE MATRIX, WITH UNIT DIAGONAL 
ELEMENTS AND NULL OTHER ELEMENTS 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A (4, 4) 

DO 5 1=1, N 
DO 5 J=1,M 
A(l , J) =0. 

IF (I.EQ.J) A(l , J) =1 . 

5 CONTINUE 
RETURN 
END 


SUBROUTINE NONLIN 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON /BLOCK1/ X(13) , Y (10) ,A(10, 10) ,Y1 (10) , IPVT(IO) ,WORK(10) , 
+ EPSI, DELTA, NC,NE,NDIM 


DO 5 1=1, NC 
CALL FUNC 

WRITE (6,10) I, (X(J) ,Y(J) , J-1,5) 
10 FORMAT (IX, '***' , 15/5 (IX, 2D15 . 7/) ) 
DO 15 J= 1 , NE 

IF (DABS (Y (J) ) .GT. EPSI) GO TO 25 
15 CONTINUE 
GO TO 105 
25 DO 35 J=1 ,NE 
35 Y1(J)=Y(J) 
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DO 45 J=1 ,NE 
DIFF=DABS (X ( J) ) *DELTA 
IF (DABS(X(J)) .LT.l.D-12) DIFF=DELTA 
XMAM=X(J) 

X(J) =X (j) -DIFF 
CALL FUNC 
X ( J) =XMAM 
DO 55 K=1,NE 

55 A(K, J) = (Y 1 (K) -Y (K) ) /DIFF 
45 CONTINUE 
DO 65 J= 1 , NE 
65 Y(J)=-Y1(J) 

C DO 205 K= 1 , NE 
C 205 WRITE (6,245) (A(K, J) , J=1 ,NE) 

C 245 FORMAT (IX, ' , 5D15 . 7) 

CALL DECOMP (NDIM.NE, A.COND, IPVT, WORK) 
CALL SOLVE (NDIM,NE,A, Y, IPVT) 

DO 75 J=1 ,NE 
75 X(J)=X(J)+Y(J) 

5 CONTINUE 
C 105 WRITE (6,20) I 
C 20 FORMAT (IX, '1=' ,12) 

105 RETURN 
END 


SUBROUTINE DECOMP (NDIM,N, A,COND, IPVT, WORK) 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A(NDIM,N) , WORK (N) , IPVT (N) 

C 

C DECOMPOSES AREAL MATRIX BY GAUSSIAN ELIMINATION, 

C AND ESTIMATES THE CONDITION OF THE MATRIX. 

C 

C -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS-, BY G. E. FORSYTHE, 
C M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

C 

C USE SUBROUTINE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEM. 

C 

C INPUT.. 

C 

C NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A 

C N ORDER OF THE MATRIX 

C A MATRIX TO BE TRIANGULARIZED 

C 

C OUTPUT.. 

C 

C A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PREMUTED 
C VERSION OF A LOWER TRIANGULAR MATRIX I~L SO THAT 

C (PERMUTATION MATRIX) *A=L*U 

C 

C COND = AN ESTIMATE OF THE CONDITION OF A. 

C FOR THE LINEAR SYSTEM A*X = B , CHANGES IN A AND B 
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C MAY CAUSE CHANGES COND TIMES AS LARGE IN X. 

C IF COND+l.O .EQ. COND , A IS SINGULAR TO WORKING 

C PRECISION. COND IS SET TO 1.0D+32 IF EXACT 

C SINGULARITY IS DETECTED. 

C 

C IPVT = THE PIVOT VECTOR 

C IPVT (K) = THE INDEX OF THE K~TH PIVOT ROW 

C IPVT (N) = (-1)** (NUMBER OF INTERCHANGES) 

C 

C WORK SPACE.. THE VECTOR WORK MUST BE DECLARED AND INCLUDED 
C IN THE CALL. ITS INPUT CONTENTS ARE IGNORED. 

C ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. 

C 

C THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY 
C DET(A) = IPVT (N) * A(1 , 1) * A(2,2) * ... * A(N,N) . 

C 

IPVT (N) =1 

IF (N. EQ. 1) GO TO 150 
NM1=N-1 

C COMPUTE THE 1-NORM OF A . 

ANORM=0 . DO 
DO 20 J=1,N 
T=O.DO 
DO 10 1=1, N 
10 T=T+DABS (A(I , J) ) 

IF (T.GT.ANORM) ANORM=T 
20 CONTINUE 

C DO GAUSSIAN ELIMINATION WITH PARTIAL 

C PIVOTING. 

DO 70 K=1 , NM1 
KP 1 =K+ 1 

FIND THE PIVOT. 

M=K 

DO 30 I=KP1,N 

IF (DABS (A(I ,K) ) .GT.DABS (A(M, K) ) ) M=I 
CONTINUE 
IPVT (K) =M 

IF (M.NE.K) IPVT (N) “-IPVT (N) 

T=A(M,K) 

A(M,K) =A(K,K) 

A (K, K) =T 

SKIP THE ELIMINATION STEP IF PIVOT IS ZERO. 
IF (T.EQ.O.DO) GO TO 70 

COMPUTE THE MULTIPLIERS. 

DO 40 I=KP1,N 
A(I,K)=-A(I,K)/T 

INTERCHANGE AND ELIMINATE BY COLUMNS. 

DO 60 J=KP1,N 
T=A(M, J) 

A(M, J)=A(K, J) 

A(K, J)=T 

IF (T.EQ.O.DO) GO TO 60 
DO 50 I=KP1,N 


C 


30 


C 


40 

C 
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50 A(I,J)-A(I,J)+A(I,K)*T 

60 CONTINUE 
70 CONTINUE 
C 

C COND = (1-NORM OF A)* (AN ESTIMATE OF THE 1-NORM OF A-INVERSE) 

C THE ESTIMATE IS OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE 
C SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS 
C OF EQUATIONS, (A- TRANSPOSE) *Y = E AND A*Z = Y WHERE E 
C IS A VECTOR OF +1 OR -1 COMPONENTS CHOSEN TO CAUSS GROWTH IN Y. 

C ESTIMATE = (1-NORM OF Z) / (1-NORM OF Y) 

C 

C SOLVE (A- TRANSPOSE) *Y = E . 

DO 100 K=1 , N 
T=0.D0 

IF (K.EQ.l) GO TO 90 
KM1=K-1 
DO 80 1=1 , KMl 
80 T=T+A(l,K)*WORK(l) 

90 EK= 1 . DO 

IF (T.LT.O.DO) EK=-1.D0 
IF (A(K,K) .EQ.O.DO) GO TO 160 
100 WORK (K) =- (EK+T) /a(K,K) 

DO 120 KB=1,NM1 
K=N-KB 
T=0.D0 
KP1=K+1 

DO 110 I=KP1,N 
110 T=T+A(I,K)*W0RK(K) 

WORK(K)=T 

M=IPVT(K) 

IF (M.EQ.K) GO TO 120 
T=WORK(M) 

WORK (M) =WORK(K) 

WORK (K)=T 
120 CONTINUE 
C 

YNORM=O.DO 
DO 130 1=1, N 

130 YNORM=YNORM+DABS (WORK (I) ) 

SOLVE A*Z = Y 

CALL SOLVE (NDIM,N, A, WORK, IPVT) 

ZNORM=O.DO 
DO 140 1=1, N 

140 ZNORM=ZNORM+DABS (WORK ( I) ) 

ESTIMATE THE CONDITION. 

COND=ANORM*ZNORM/ YNORM 
IF (COND. LT. 1 .DO) COND=1.DO 
RETURN 

C 1-BY-l CASE.. 

150 COND= 1 . DO 

IF (A(l,l) .NE.O.DO) RETURN 
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EXACT SINGULARITY 

C0ND=1 . 0D32 

RETURN 

END 

SUBROUTINE SOLVE (NDIM, N , A, B , IPVT) 


IMPLICIT REAL*8(A-H,0-Z) 
DIMENSION A (NDIM, N) ,B (N) , IPVT (N) 


SOLVES A LINEAR SYSTEM, A*X = B 

DO NOT SOLVE THE SYSTEM IF DECOMP HAS DETECTED SINGULARITY. 


-COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS", BY G. E. FORSYTHE, 
M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 


INPUT. . 

NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A 
N = ORDER OF MATRIX 

A = TR I ANGULAR IZED MATRIX OBTAINED FROM SUBROUTINE DECOMP 

B = RIGHT HAND SIDE VECTOR 

IPVT = PIVOT VECTOR OBTAINED FROM DECOMP 


OUTPUT . . 

B = SOLUTION VECTOR, X 


DO THE FORWARD ELIMINATION. 

IF (N.EQ.l) GO TO 50 
NM1=N-1 
DO 20 K=1 ,NM1 
KPl=K+l 
M=IPVT(K) 

T=B (M) 

B (M) =B (K) 

B (K) =T 

DO 10 I=KP1,N 
10 B(l)=B(l)+A(l,K)*T 

20 CONTINUE 

NOW DO THE BACK SUBSTITUTION. 

DO 40 KB=1,NM1 
KM1=N-KB 
K=KM1+1 

B(K) =B (K) / A (K , K) 

T=-B (K) 

DO 30 1=1, KM1 
30 B(l)=B(l)+A(l,K)*T 

40 CONTINUE 
50 B(1)=B(1)/A(1,1) 

RETURN 

END 
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(2 _ _ ******************************************************* 

c 

C PURPOSE 
C 

C THIS PROGRAM IS USED TO CALCULATE THE TRANSMISSION ERRORS OF A 
C HELICAL PINION AND A HELICAL GEAR IN MESHING WHEN THEIR AXES ARE 
C DEFORMED BY INTERACTING FORCE (BOTH PINION AND GEAR ARE NOT CROWNED) 
C 

C NOTE 
C 

C THIS PROGRAM IS WRITTEN IN FORTRAN 77. IT CAN BE COMPILED BY V 
C COMPILER IN IBM MAINFRAME OR FORTRAN COMPILER IN VAX SYSTEM. 

C 


IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,MUO,KSIC 

DIMENS ION Z (99) , ANGLE (99) , ERROR (99) , ERR (99) , W (99) , RPR (99) , RGR (99) , 
+ Sl(4,4),S2(4,4),ZP(99) 

COMMON /BLOCK1/ X(12) , Y(10) ,A(10, 10) , Y1 (10) , IPVT(IO) ,WORK(10) , 

+ EPSI, DELTA, NC,NE,NDIM,NCTL,CX2 

COMMON /BLOCK2/ S (4 , 4) , C, RP ,RG, CK, SK, CB , SB , CKC , KSIC , ZG, TB , ZNPF , 

+ C0EG1 , C0EG2,CB1 , SB1 , TB1 

C 

C DEFINE PARAMETERS USED BY PROGRAMS 
C 

C (1) IN ANG LP ARE UNIT NUMBERS ASSIGNED TO THE INPUT AND OUTPUT 
C DEVICES 


IN=5 


LP=6 

C (2) NDBUG IS USE TO CONTROL THE AUXILIARY OUTPUT FOR DEBUGGING 
NDBUG= 1 

C (3) NC IS THE UPPER LIMITATION OF REPEATATION FOR SOLVING NONLINEAR 
C EQUTIONS; 

C EPSI IS THE CLEARANCE OF FUNCTION VALUES WHEN THE FUNCTIONS 

C IS CONSIDERED AS SOLVED (ALL FUNTIONS HAVE FORMS OF F(X)=0); 

C DELTA IS THE RELATIVE DIFFERENCE FOR TAKING DERIVATIVES 

C NC, EPSI AND DELTA MAY BE CHANGED WHEN SOLUTIONS ARE DIVERGENT 

C OR LESS ACCURATE 

NC=100 


DELTA=l.D-4 


EPSI=1.D-13 

C (4) OTHER PARAMETERS (DON’T CHANGE) 
NDIM=10 

DR-DATAN ( 1 . DO) /45 . DO 
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C DEFINE INPUT PARAMTERS OF PROBLEM (USE INCH AS UNIT OF LENGTH) 

C (1) PINION AND GEAR: PN=DIAMETRAL PITCH; NP=PINI0N TOOTH NUMBER; 
C RMPG- TOOTH NUMBER RATIO (GEAR TOOTH NO./NP); 

C KSIN-PRESSURE ANGLE IN NORMAL SECTION; 

C BETAP=HELIX ANGLE OF PINION AND GEAR; 

C PAD=HEIGHT OF ADDENDUM OF PINION; 

C GAD-HEIGHT OF ADDENDUM OF GEAR; 

C ZG-GEAR TOOTH LENGTH (PINION TOOTH LENGTH IS LONGER) 

C COE=COEFF. OF CENTRAL DISTANCE (USUALLY COE-1.) 

PN-10.D0 
NP-20 
RMPG-2.D0 
KS IN-20. D0*DR 
BETAP-15 . D0*DR 
PAD=1 . 0/PN 
GAD=1 . 0/PN 
ZG= 5./PN 
COE=l . 000D0 

C (2) SHAFT DEFORMATION: 

C NSIM-MODEL ID NO. (1-SIMPLIFIED DEFORMATION MATRIX; 

C 2-UNSIMPLIFIED DEFORMATION MATRIX) 

C RLAMP-PINION SHAFT ROTATION 

C RLAMG-GEAR SHAFT ROTATION 

C RVP=PINION SHAFT DEFLECTION 

C RVG-GEAR SHAFT DEFLECTION 

C 

NSIM=1 

RLAMP= 2./60.*DR 
RLAMG- 2./60.*DR 
RVP=0. 0125 
RVG-0.0125 

C (3) OUTPUT: FEE 1= INCREMENT OF ROTATION ANGLE OF PINION (DEGREEE) 
FEEI=l.OD0*DR 
C 

C DESCRIPTION OF OUTPUT PARAMERTERS 
C 

C FEE 1 -ROTATION ANGLE OF PINION 

C FEE2-ROTATION ANGLE OF GEAR 

C RP-RADIUS OF PINION CONTACT POINT 
C RG-RADIUS OF GEAR CONTACT POINT 
C 

C FIND AUXILIARY VALUES FOR CALCULATION 

C 

C 

C DEFINE USEFUL CONSTANTS AND PARAMETERS FOR PINION AND GEAR 
DELTAB-DR/60 . * ( 0.D0) 

BETAP 1-BETAP+DELTAB 
CLP-DCOS (RLAMP) 

SLP-DSIN (RLAMP) 

CLG-DCOS (RLAMG) 

SLG-DSIN (RLAMG) 

WRITE (LP, 4) CLP , SLP , CLG , SLG 
4 FORMAT (IX, 4F15 . 7) 

SK-DSIN (KSIN) 
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CK=DCOS (KSIN) 

SB=DSIN (BETAP) 

CB=DCOS (BETAP) 

TB=SB/CB 
SB1=DSIN (BETAP1) 

CBl=DCOS (BETAP 1) 

TB1=SB1/CB1 
KSIC=DATAN (SK/CK/CB) 

CKC=DCOS(KSIC) 

SKC=DSIN (KSIC) 

PT=PN*CB 
RP=NP/2./PT 
RPA=RP+PAD 
RPATOL=RPA-0 . 0005D0 
RG=RP*RMPG 
RGA=RG+GAD 
RGATOL=RGA-0 . 0005DO 
WRITE (LP , 56) RPATOL.RGATOL 
56 FORMAT (IX, ' &&&&&& ' , 'RPA= ' , F15 . 7, 5X, 'RGA= 1 , F15 . 7) 

COEGl=l . +SK*SK*TB*TB 
COEG2=RG/COEGl 
ZNPF=CK*SB1 
C= (RP+RG) *COE 
CALL INTMAT(S, A,4) 

CALL INTMAT (S 1 , 4 , 4) 

CALL INTMAT (S2, 4, 4) 

DO 505 LL=1 , 2 
NSIM=LL 

IF (NSIM.EQ.l) THEN 
WRITE (LP,500) 

500 FORMAT (1H1,///,1X, 

+ 'THE CASE OF SIMPLIFIED DEFORMATION MATRIX OF GEAR AXES') 

S ( 1 , 3) = (RLAMP+RLAMG) *CKC 
S (2 , 3) = (RLAMP+RLAMG) *SKC 
S (3, 1) =~S (1 , 3) 

S(3,2)— S(2,3) 

S (1 , 4) = (RVP+RVG) *CKC 
S (2 , 4) - (RVP+RVG) *SKC+C 
S (3 , 4) =~C*RLAMP*SKC 
ELSE 

WRITE (LP , 501) 

501 FORMAT (1H1,///,1X, 

+ 'THE CASE OF UNSIMPLIFIED DEFORMATION MATRIX OF GEAR AXES') 

S 1 ( 1 , 1 ) = 1 . +CKC*CKC* (CLP- 1 . ) 

SI (1 , 2) =SKC*CKC* (CLP-1 . ) 

SI (1,3)=CKC*SLP 
SI (2 , 1) =S1 (1 , 2) 

SI (2 , 2) =1 . +SKC*SKC* (CLP-1 . ) 

SI (2 , 3) =SKC*SLP 
SI (3,1)— SI (1,3) 

SI (3, 2) —SI (2,3) 

SI (3,3)=CLP 
SI (1,4) =RVP*CKC*CLP 
SI (2, 4) =RVP*SKC*CLP 
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S1(3,4)=-RVP*SLP 

S2 (1 , 1) =1 . +CKC*CKC* (CLG-1 . ) 

S2 (1 , 2) =SKC*CKC* (CLG-1 . ) 

S2 ( 1 , 3) =CKC*SLG 
S2(2, 1) =S 1 (1,2) 

S2 (2 , 2) =1 . +SKC*SKC* (CLG- 1 . ) 

S2 (2 , 3) =SKC*SLG 
S2(3, 1 ) =— S 1 (1,3) 

S2 (3 , 2) =~S1 (2 , 3) 

S2 (3,3) =CLG 
S2 (1 , 4) =RVG*CKC 
S2 (2 , 4) =RVG*SKC+C 
S2 (3 , 4) =0. 

DO 12 MMI=1 , 3 
DO 12 MMJ=1,4 
S (MMI ,MMJ) =0. 

DO 12 MMK=1,4 

S (MMI , MM J) =S (MMI , MM J) +S 1 (MMI , MMK) *S2 (MMK,MMJ) 
12 CONTINUE 
END IF 

NCOEF=IDINT (360. *DR/ FEE I /FLOAT (NP) +0 . 5) 

N=3*NCOEF 

NHALF= (N+l)/2 

XIN=0 . 

DO 205 L=1 , 2 
LSGN= (-1) **l 
NCTL=0 
NE=4 

DO 5 1=1, NE 
5 X(l)=0.D0 
X(10)=0. 

X(11)=0. 

DO 15 1=1 , NHALF 

IF (NCTL.EQ.2) GOTO 205 

LI=NHALF+LSGN*(I-1) 

IF (L.EQ.l) NMIN=LI 
IF (L.EQ.2) NMAX=LI 
X (7) = LSGN*FEE I *FLOAT (1-1) 

IF (X ( 1 0) . LE . RPATOL) THEN 
X(4) =-X(7) 

ELSE 

NE=3 

NCTL=NCTL+1 
END IF 

IF (X (11) .LT.RGATOL) THEN 
X (2) =X (7) /RMPG 
ELSE 

NCTL=NCTL+1 

CX2=X(2) 

END IF 

X (3) = (ZG-RP*X (4) *SK*SK*TB) /CB/COEG1 

CALL NONLIN 

X(8)=X(1)+X(2) 

C FIND INITIAL VALUE OF X(8) 
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IF (L.EQ.l.AND.I.EQ.l) THEN 
XIN=X(8) 

GO TO 15 
END IF 

X (8) =X (8) -XIN 
W(LI)=X(7)/DR 
Z(LI)=X(8)/DR 

ERR (LI) = (X (8) -X (7) /RMPG) *3600 . DO/DR 
RPR (LI) =X (10) 

RGR (LI) =X (1 1) 

ZP(LI)=X(12) 

WRITE (LP ,21) LI,W(LI) ,Z(LI) , ERR (LI) ,RPR(Ll) ,RGR(Ll) 

21 FORMAT (IX, 12 , IX, 5F15 . 7 , F15 . 7) 

15 CONTINUE 
205 CONTINUE 

WRITE (LP , 10) 

10 FORMAT (////8X, 'FEEl(D) ’ ,8X, ’FEE2(D) ' ,8X, 'K-ERROR(S) ' ,5X, 

+ 'RP ' , 13X, ' RG' , F15. 7/) 

DO 25 I=NMIN,NMAX 

WRITE (LP , 20) W(I) ,Z(I),ERR(I) ,RPR(I),RGR(I) 

20 FORMAT (1X,5F15 . 7, F15 . 7) 

25 CONTINUE 

NT=NMAX-NCOEF 

FII=FEEl/DR/2 . D0*FLOAT (NCOEF) 

WRITE (LP ,80) 

80 FORMAT (//, ' FIND THE WORKING RANGE FOR ONE TOOTH: ' , F15 . 7/) 
DO 55 I=NMIN,NT 

X (7) =FEEI*FLOAT (I- (N+l) /2) /DR/2 . DO 

X(8) =X(7)+FII 

KK=I+NCOEF 

ANGLE (i) = Z (KK) -Z (I) 

ERROR (I) = (ANGLE (i) -FII) *3600. DO 

WRITE (LP ,60) X (7) ,X (8) .ANGLE (I) .ERROR (I) 

60 FORMAT (IX, ' (' ,F7.2, ' ' ,F7.2, ') : ' , F15 . 7, F15 . 7) 

55 CONTINUE 

DO 95 I=NMIN,NT 
ATEMP2=ERROR (I) 

IF (I.NE.NMIN) THEN 
IF (ATEMPl*ATEMP2.LE.0.D0) GOTO 105 
END IF 

ATEMP1=ATEMP2 
95 CONTINUE 

WRITE (LP , 160) 

160 FORMAT (//IX, 'MESHING IS DISCONTINUOUS') 

GO TO 505 

105 IF (DABS (ATEMP1) .LT. DABS (ATEMP2)) 1-1*1 
EMAX=0. 

EMIN=0. 

NTEMP“NCOEF+ 1 
DO 135 J=1 , NTEMP 
KS=I+J~1 
ET=ERR (KS) 

IF (ET. LT. EMIN) EMIN=ET 
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IF (ET.GT.EMAX) EMAX=ET 
135 CONTINUE 

ET=EMAX-EMIN 

KK=I+NCOEF 

WRITE (LP ,170) Z(I) ,2(KK) ,ET 

170 FORMAT (//IX , 'WORKING RANGE FOR ONE TOOTH: ',¥1.2,' ' ,F7.2/ 

+ IX, 'THE MAXIMUM KINEMATIC ERROR: ',Fl5.7,' (S)',I2) 

505 CONTINUE 
STOP 
END 
C 

SUBROUTINE FUNC 
C 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION KSIN,MU,MUO,KSIC 

COMMON /BLOCK1/ X (12) , Y (10) , A (10, 10) , Yl (10) , IPVT (10) , WORK (10) , 

+ EPSI , DELTA, NC , NE, NDIM, NCTL, CX2 

COMMON /BLOCK2/ S (4 , 4) , C ,RP , RG, CK, SK, CB , SB , CKC.KSIC , ZG, TB , ZNPF, 
+ COEG1 ,COEG2, CB1 , SB1 , TBl 

C X(4)=FEEP; X (7) =FEE1 
X (6) =X(4)+X(7) 

CF1FP=DC0S(X (6)) 

SF1FP=DSIN (X (6) ) 

AP=X (3) *SB 1+RP*X (4) 

XPF=-CK*CK/CKC*DCOS (X(6) -KSIC) *AP+RP*SF1FP 
YPF=CK*CK/CKC*DSIN (X (6) -KSIC) *AP+RP*CF1FP 
ZPF=X(3) *CB1+AP*SK*SK*TB1 
XNPF=CK*CB1*CF1FP+SK*SF1FP 
YNPF=-CK*CB1*SF1FP+SK*CF1FP 
C ZNPF=CK*SB1 

CF2FG=DCOS (X(l)) 

SF2FG=DSIN (X (1) ) 

CF2FGK-DCOS (X(1)+KSIC) 

SF2FGK=DSIN (X (1) +KSIC) 

AG1= (-ZG*TB+RG*X (2) ) /COEG1 
RFG1=CK*CK/CKC*CF2FGK*AG1+RG*SF2FG 
RFG2=CK*CK/CKC*SF2FGK*AG1-RG*CF2FG 
C RFG3=ZG 

RTFG1=CK*CK/CKC* ( AG1*SF2FGK+C0EG2*CF2FGK) -RG*CF2FG 
RTFG2=CK*CK/CKC* (-AG1*CF2FGK+C0EG2*SF2FGK) -RG*SF2FG 
C RTFG3=0. 

XGF=RFG1*S (1,1) +RFG2*S (1,2) +ZG*S ( 1 , 3) +S (1 , A) 

YGF=RFG1*S (2 , 1) +RFG2*S (2 , 2) +ZG*S (2 , 3) +S (2 , A) 

ZGF=RFG1*S (3 , 1) +RFG2*S (3 , 2) +ZG*S (3 , 3) +S (3 , 4) 

XTGF=RTFG1*S (1 , 1) +RTFG2*S (1 , 2) 

YTGF=RTFG1*S (2,1) +RTFG2*S (2 , 2) 

ZTGF=RTFG1*S (3 , 1) +RTFG2*S (3 , 2) 

C WRITE (6,100) RTFG1 ,RTFG2, AG1 ,C0EG2, COEG1 ,CF2FGK, SF2FGK 

C WRITE (6,100) X(l) ,X(2) ,X(3) ,X(4) ,X(5) 

C WRITE (6,100) XPF , YPF, ZPF.XGF, YGF , ZGF 

C 100 FORMAT (IX, "LVLVL' .8E15.7) 

C WRITE (6,100) XNPF , YNPF , ZNPF , XTGF , YTGF , ZTGF 

Y(1)=XPF-XGF 
Y(2)=YPF-YGF 
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Y(3)=ZPF-ZGF 
IF (NCTL.NE.O) THEN 
Y (A) =X (2) -CX2 
ELSE 

Y (4) =XNPF*XTGF+YNPF*YTGF+ZNPF*ZTGF 
END IF 

X (10) =DSQRT (XPF*XPF+YPF*YPF) 

X (1 1) =DSQRT (RFG1*RFG1+RFG2*RFG2) 

X (12) =ZPF 

WRITE (6,20) (Y(II) ,11=1,4) 

20 FORMAT (IX, ' $$$$ ' ,6F15. 7) 

RETURN 

END 

SUBROUTINE INTMAT (A,N,M) 

THIS SUBROUTINE IS USED TO INITIATE THE MATRIX, WITH UNIT DIAGONAL 
ELEMENTS AND NULL OTHER ELEMENTS 
IMPLICIT REAL*8(A~H,0-Z) 

DIMENSION A (4, 4) 

DO 5 1=1, N 
DO 5 J=1,M 
A(I, J)=0. 

IF (I.EQ.J) A(I,J)=1. 

5 CONTINUE 
RETURN 
END 


SUBROUTINE NONLIN 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON /BLOCK1/ X (12) , Y (10) , A(10 , 10) , Yl (10) , IPVT (10) , WORK (10) , 
+ EPS I , DELTA , NC , NE , NDIM , NCTL , CX2 

DO 5 1=1, NC 
CALL FUNC 

WRITE (6,10) I, (X(j) , Y ( J) , J=1 , 4) 

10 FORMAT (IX, '***' , 15/5 (IX, 2D15 . 7/) ) 

DO 15 J=1 , NE 

IF (DABS (Y ( J) ) . GT. EPSI) GO TO 25 
15 CONTINUE 
GO TO 105 
25 DO 35 J=1,NE 
35 Yl ( J) =Y ( J) 

DO 45 J=1 ,NE 
DIFF=DABS (X (J) ) *DELTA 
IF (DABS(X(J)) .LT.l.D-12) DIFF=DELTA 
XMAM=X ( J) 

X(J)=X(J)-DIFF 
CALL FUNC 
X(J) =XMAM 
DO 55 K= 1 , NE 

55 A(K,J) = (Y1(K)-Y(K))/DIFF 
45 CONTINUE 
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DO 65 J-l.NE 
65 Y(J)=-Y1(J) 

C DO 85 K=1 ,NE 

C 85 WRITE (6,104) (A(K, J) , J-l.NE) ,Y(K) 

C 104 F0RMAT(1X, 'A' ,2X,5E15.7) 

CALL DECOMP (NDIM,NE,A,COND, I PVT, WORK) 
CALL SOLVE (NDIM.NE, A, Y, IPVT) 

DO 75 J= 1 , NE 
75 X(J)-X(J)+Y(J) 

5 CONTINUE 
105 RETURN 
END 


SUBROUTINE DECOMP (NDIM.N, A, COND, IPVT, WORK) 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A (NDIM.N) .WORK (N) .IPVT (N) 

C 

C DECOMPOSES AREAL MATRIX BY GAUSSIAN ELIMINATION, 

C AND ESTIMATES THE CONDITION OF THE MATRIX. 

C 

C -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS-, BY G. E. FORSYTHE, 
C M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

C 

C USE SUBROUTINE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEM. 

C 

C INPUT.. 

C 

C NDIM - DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A 

C N ORDER OF THE MATRIX 

C A MATRIX TO BE TRIANGULARIZED 

C 

C OUTPUT.. 

C 

C A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PREMUTED 
C VERSION OF A LOWER TRIANGULAR MATRIX I-L SO THAT 

C (PERMUTATION MATRIX) *A=L*U 

C 

C COND = AN ESTIMATE OF THE CONDITION OF A. 

C FOR THE LINEAR SYSTEM A*X = B , CHANGES IN A AND B 

C MAY CAUSE CHANGES COND TIMES AS LARGE IN X. 

C IF COND+l.O .EQ. COND , A IS SINGULAR TO WORKING 

C PRECISION. COND IS SET TO 1.0D+32 IF EXACT 

C SINGULARITY IS DETECTED. 

C 

C IPVT = THE PIVOT VECTOR 

C IPVT (K) = THE INDEX OF THE K-TH PIVOT ROW 

C IPVT (N) = (-1)** (NUMBER OF INTERCHANGES) 

C 

C WORK SPACE.. THE VECTOR WORK MUST BE DECLARED AND INCLUDED 
C IN THE CALL. ITS INPUT CONTENTS ARE IGNORED. 

C ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. 
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THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY 
DET(A) = IPVT(N) * A(l,l) * A(2,2) * ... * A(N,N) . 

IPVT (N) =1 

IF (N.EQ.l) GO TO 150 
NM1-N-1 

COMPUTE THE 1-NORM OF A . 

ANORM=0 . DO 
DO 20 J=1,N 
T=0.D0 
DO 10 1=1, N 
10 T=T+DABS (A (I , J) ) 

IF (T.GT. ANORM) ANORM=T 
20 CONTINUE 

DO GAUSSIAN ELIMINATION WITH PARTIAL 
PIVOTING. 

DO 70 K=1 ,NM1 
KP1=K+1 

FIND THE PIVOT. 

M=K 

DO 30 I=KP1,N 

IF (DABS (A(l ,K) ) .GT.DABS (A(M,K) ) ) M=I 
30 CONTINUE 
IPVT (K) =M 

IF (M.NE.K) IPVT(N)=-IPVT(N) 

T=A(M,K) 

A (M , K) =A(K,K) 

A(K,K)=T 

SKIP THE ELIMINATION STEP IF PIVOT IS ZERO. 
IF (T.EQ.O.DO) GO TO 70 

COMPUTE THE MULTIPLIERS. 

DO 40 I=KP1,N 
40 A(I,K)=-A(I,K)/T 

C INTERCHANGE AND ELIMINATE BY COLUMNS. 

DO 60 J=KP 1 , N 
T=A(M, J) 

A(M, J) =A(K, J) 

a(k, j) =t 

IF (T.EQ.O.DO) GO TO 60 
DO 50 I=KP1,N 

50 A(I, J)=A(I, J)+A(I,K)*T 
60 CONTINUE 
70 CONTINUE 
C 

C COND = (1-NORM OF A)* (AN ESTIMATE OF THE 1-NORM OF A-INVERSE) 

C THE ESTIMATE IS OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE 
C SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS 
C OF EQUATIONS, (A- TRANSPOSE) *Y = E AND A*Z = Y WHERE E 
C IS A VECTOR OF +1 OR -1 COMPONENTS CHOSEN TO CAUSS GROWTH IN Y. 

C ESTIMATE = (1-NORM OF Z) / (1-NORM OF Y) 

C 

C SOLVE (A- TRANSPOSE) *Y = E . 
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DO 100 K=1,N 
T=0.D0 

IF (K.EQ.l) GO TO 90 
KM1=K-1 
DO 80 1=1, KM1 
80 T=T+A(l,K)*WORK(l) 

90 EK= 1 . DO 

IF (T.LT.0.D0) EK=- 1 . DO 
IF (A(K,K) .EQ.0.D0) GO TO 160 
100 WORK (K) =~ (EK+T) /A (K, K) 

DO 120 KB=1,NM1 
K=N-KB 
T=0.D0 
KP1=K+1 

DO 110 I=KP1 , N 

no t=t+a(i,k)*work(k) 

WORK (K) =T 
M=IPVT (K) 

IF (M.EQ.K) GO TO 120 
T=WORK (M) 

WORK (M) =WORK(K) 

WORK (K) =T 
120 CONTINUE 
C 

YNORM=O.DO 
DO 130 1=1, N 

130 YNORM= YNORM+DAB S (WORK (I) ) 

C 

C SOLVE A*Z = Y 

CALL SOLVE (NDIM.N, A, WORK, IPVT) 

C 

ZN0RM=0.D0 
DO 140 1=1, N 

140 ZNORM=ZNORM+DAB S (WORK ( I ) ) 


ESTIMATE THE CONDITION. 

COND= ANORM*ZNORM/ YNORM 
IF (COND.LT. 1.D0) COND-1.D0 
RETURN 

1-BY-l CASE.. 


150 COND= 1 . DO 

IF (A(l,l) .NE.O.DO) RETURN 


EXACT SINGULARITY 

160 COND=l . 0D32 
RETURN 
END 

SUBROUTINE SOLVE (NDIM.N, A, B, IPVT) 

C 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A(NDIM,N) ,B (N) , IPVT (N) 

C 

C SOLVES A LINEAR SYSTEM, A*X = B 

C DO NOT SOLVE THE SYSTEM IF DECOMP HAS DETECTED SINGULARITY. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 


-COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS-, BY G. E. FORSYTHE, 
M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 


INPUT. . 

NDIM - DECLARED ROW DIMENSION OF ARRAY CONTAINING A 
N = ORDER OF MATRIX 

A = TRIANGULAR IZED MATRIX OBTAINED FROM SUBROUTINE DECOMP 

B = RIGHT HAND SIDE VECTOR 

IPVT = PIVOT VECTOR OBTAINED FROM DECOMP 

OUTPUT. . 

B = SOLUTION VECTOR, X 


DO THE FORWARD ELIMINATION. 

IF (N.EQ.l) GO TO 50 
NM1=N-1 
DO 20 K=1 ,NM1 
KP1=K+1 
M=IPVT (K) 

T=B (M) 

B (M) =B (K) 

B(K)=T 

DO 10 I=KP1,N 
10 B (I) =B (i) +A(l ,K) *T 

20 CONTINUE 

NOW DO THE BACK SUBSTITUTION. 

DO 40 KB-l.NMl 
KM1=N-KB 
K-KM1+1 

B(K)=B (K) /A(K,K) 

T=-B (K) 

DO 30 1=1, KM1 
30 B(I)=B(I)+A(I,K)*T 

40 CONTINUE 
50 B(1)=B(1)/A(1,1) 

RETURN 

END 
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